はじめに
前回、pythonを用いてFFT解析を行う手法について解説しました。
今回は、フィルタ処理のやり方を見ていきましょう。
FFT解析ができれば、ノイズである周波数を取り除いたり、エイリアシングを防ぐために高周波成分を取り除いたりするようなフィルター処理が必要になってきます。
さっそくやっていきましょう。
ローパスフィルタ
今回は題材として低周波成分のみ抽出するローパスフィルタを作成してみます。
昔MATLABを使用していましたが、その際によく使っていたのがバターワースフィルタというもので、scipyにもMATLABチックに実装されてますのでそれを使用していきます。
バターワースフィルタは、その次数に伴って挙動が変わります。
下記に次数の違うバターワースローパスフィルタの周波数応答を示します。
カットオフ周波数は全部同じ100Hzで、次数を1~4まで変えて応答を見ています。
これを見るに、すべてのフィルタは100Hzの手前からゲインが下がりはじめ、周波数が大きくなるほど0に漸近するような応答をしています。
次数による違いは100Hz前後での変化の傾きが異なります。
次数が高いほうが急峻に変化するので、精度の高いフィルタとなることがわかります。
バターワースフィルタの実装
scipy.signalに関連のパッケージが実装されておりますのでそれを使用します。
scipy.signal.butter()でバターワースフィルタの設計。
scipy.signal.filtfilt()でデータにフィルタをかける処理を行います。
バターワースフィルタの周波数応答
先ほどの図を描いたコードは以下です。
カットオフ周波数はナイキスト周波数(サンプリング周波数の半分)で割って正規化した値をscipy.signal.butterの引数に指定してあげる必要があります。
import numpy as np import matplotlib.pyplot as plt from scipy import signal def butter_lowpass(lowcut, fs, order=4): '''バターワースローパスフィルタを設計する関数 ''' nyq = 0.5 * fs low = lowcut / nyq b, a = signal.butter(order, low, btype='low') return b, a # Plotting the frequency response. fs = 2000 # sampling frequency cutoff = 100 # cutoff frequency plt.figure() # 次数1~4まで繰り返し描画 for i in range(1, 5): b, a = butter_lowpass(cutoff, fs, order=i) w, h = signal.freqz(b, a) plt.plot(0.5*fs*w/np.pi, np.abs(h), label='order='+str(i)) plt.axvline(cutoff, color='k', label='cutoff={}Hz'.format(cutoff)) plt.xlim(0, 0.5*fs) plt.title("Lowpass Filter Frequency Response") plt.xlabel('Frequency [Hz]') plt.legend(loc='best') plt.grid()
周波数応答を求めるのはscipy.signal.freqz()パッケージを使用すれば一発です。
freqzの返り値は周波数wとその時の応答hです。
wの値はナイキスト周波数で正規化されているのでナイキスト周波数をかけてあげる必要があります。
ローパスフィルターをかけてみる
実際にフィルターをかけてみましょう。
サンプルデータの準備
サンプルとして以下データを用意します。
import numpy as np import matplotlib.pyplot as plt from scipy import signal N = 1024 # サンプル数 dt = 0.0005 # サンプリング周期 [s] fs = 1 / dt # サンプリング周波数 [Hz] f1, f2, f3 = 30, 432, 604 # サンプルデータの周波数 [Hz] t = np.arange(0, N * dt, dt) # 時間 [s] x = 1.5 * np.sin(2 * np.pi * f1 * t) \ + 0.5 * np.sin(2 * np.pi * f2 * t) \ + 0.7 * np.sin(2 * np.pi * f3 * t) # データ fig, ax = plt.subplots() ax.plot(t, x) ax.set_xlabel("Time [s]") ax.set_ylabel("Signal") ax.set_xlim([0, 0.1]) ax.grid() plt.show()
周波数30Hzのデータに高周波のノイズ(432,604Hz)が乗っているというサンプルです。
以前やった方法でFFT解析を行います。
FFT解析の方法は下記記事を参照ください。
# 高速フーリエ変換 def calc_amp(data, fs): '''フーリエ変換して振幅スペクトルを計算する関数 ''' N = len(data) window = signal.hann(N) F = np.fft.fft(data * window) freq = np.fft.fftfreq(N, d=1/fs) # 周波数スケール F = F / (N / 2) # フーリエ変換の結果を正規化 F = F * (N / sum(window)) # 窓関数による振幅減少を補正する Amp = np.abs(F) # 振幅スペクトル return Amp, freq Amp, freq = calc_amp(x, fs) fig, ax = plt.subplots() ax.plot(freq[:N//2], Amp[:N//2]) ax.set_xlabel("Frequency [Hz]") ax.set_ylabel("Amplitude") ax.grid() plt.show()
思惑通り30,432,604Hzの所にピークがたっており、サンプルデータはこれらの周波数成分のたしあわせであることがわかります。
ここで、ローパスフィルタを用いてノイズである100Hz以上の周波数成分をフィルタリングしてみます。
フィルタ用の関数定義
フィルタを適用するのに以下2つの関数を作っておきます。
ナイキスト周波数で正規化したりいろいろ面倒なので関数化してすぐ使えるようにしておくと楽です。
from scipy import signal def butter_lowpass(lowcut, fs, order=4): '''バターワースローパスフィルタを設計する関数 ''' nyq = 0.5 * fs low = lowcut / nyq b, a = signal.butter(order, low, btype='low') return b, a def butter_lowpass_filter(x, lowcut, fs, order=4): '''データにローパスフィルタをかける関数 ''' b, a = butter_lowpass(lowcut, fs, order=order) y = signal.filtfilt(b, a, x) return y
フィルタリング
先程作った関数を用いれば簡単に実装できます。
生データにフィルタを適応した結果がこちらです。
ソースコードは以下。
import numpy as np import matplotlib.pyplot as plt from scipy import signal def butter_lowpass(lowcut, fs, order=4): '''バターワースローパスフィルタを設計する関数 ''' nyq = 0.5 * fs low = lowcut / nyq b, a = signal.butter(order, low, btype='low') return b, a def butter_lowpass_filter(x, lowcut, fs, order=4): '''データにローパスフィルタをかける関数 ''' b, a = butter_lowpass(lowcut, fs, order=order) y = signal.filtfilt(b, a, x) return y N = 1024 # サンプル数 dt = 0.0005 # サンプリング周期 [s] fs = 1 / dt f1, f2, f3 = 30, 432, 604 # 周波数 [Hz] t = np.arange(0, N * dt, dt) # 時間 [s] x = 1.5 * np.sin(2 * np.pi * f1 * t) \ + 0.5 * np.sin(2 * np.pi * f2 * t) \ + 0.7 * np.sin(2 * np.pi * f3 * t) # データ y = butter_lowpass_filter(x, 100, fs, order=4) fig, ax = plt.subplots() ax.plot(t, x, label='raw data') ax.plot(t, y, label='filtered') ax.set_xlabel("Time [s]") ax.set_ylabel("Signal") ax.set_xlim([0, 0.1]) ax.grid() plt.legend(loc='best') plt.show()
ついでにフィルタ前後でのFFTの結果を比較すると下記のようになります。
きちんと高周波成分が除外されて欲しい低周波のデータのみ抜き出すことが出来ました。
ソースコードは以下。
# 高速フーリエ変換 def calc_amp(data, fs): '''振幅スペクトルを計算する関数 ''' N = len(data) window = signal.hann(N) F = np.fft.fft(data * window) freq = np.fft.fftfreq(N, d=1/fs) # 周波数スケール F = F / (N / 2) # フーリエ変換の結果を正規化 F = F * (N / sum(window)) # 窓関数による振幅減少を補正する Amp = np.abs(F) # 振幅スペクトル return Amp, freq Amp, freq = calc_amp(x, fs) fig, ax = plt.subplots(1,2) ax[0].plot(freq[:N//2], Amp[:N//2]) ax[0].set_xlabel("Frequency [Hz]") ax[0].set_ylabel("Amplitude") ax[0].set_title('raw data') ax[0].grid() Amp_filt, freq_filt = calc_amp(y, fs) ax[1].plot(freq_filt[:N//2], Amp_filt[:N//2]) ax[1].set_xlabel("Frequency [Hz]") ax[1].set_ylabel("Amplitude") ax[1].set_title('filtered') ax[1].grid() plt.tight_layout() plt.show()
