Pythonによるデータ処理8 ～ フィルタ処理5
今回は前回紹介したバタワース、チェビシェフ、楕円特性フィルタを信号データに適用した際の周波数数成分の変化について紹介する。
1. 準備
設計したフィルタを信号データに適用する方法は以前紹介したが、今回はFFTを用いてフィルタ適用前後の周波数成分がどのように変化するか調べる。
まず、バタワースでバンドパスフィルタを設計して信号データに適用する
import numpy as np from scipy import signal from matplotlib import pyplot as plt samplerate=2500 fp = np.array([499.5,500.5]) # 通過域端周波数[Hz] fs = np.array([499,501]) # 阻止域端周波数[Hz] gpass = 3 # 通過域端最大損失[dB] gstop = 30 # 阻止域端最小損失[dB] fn = samplerate / 2 #ナイキスト周波数 wp = fp / fn #ナイキスト周波数で通過域端周波数を正規化 ws = fs / fn N, Wn = signal.buttord(fp, fs, gpass, gstop,True) #オーダーとバターワースの正規化周波数を計算 b, a = signal.butter(N, Wn, "band",True) #フィルタ伝達関数の分子と分母を計算 w, h1 = signal.freqs(b, a, np.logspace(1, 3, 100000)) fig, ax1 = plt.subplots() ax1.set_title('Butterworth response') ax1.plot(w, 20 * np.log10(abs(h1)), 'b') ax1.set_ylabel('Amplitude [dB]', color='b') ax1.set_xlabel('Frequency [rad/sample]') ax1.set_ylim([-150,10]) ax2 = ax1.twinx() angles = np.unwrap(np.angle(h1)) ax2.plot(w, angles*180/np.pi, 'g') ax2.set_ylabel('Angle (Degrees)', color='g') ax2.grid() plt.xlim([490,510])
今回は低域側は499.5Hzで通過、499Hzで阻止、高域側は500.5Hzまで通過、501Hzで阻止する狭帯域のバンドパスフィルタとする。 阻止域端最小損失は30dBとしたが、これ以上大きくするとフィルタの応答曲線がいびつになる。バタワースバンドフィルタの応答は以下のようになる。
信号データはサンプリングレート2500で10秒間のガウシアンランダムノイズを使用する。
x = np.arange(0, 25000) / samplerate data = np.random.normal(loc=0, scale=1, size=len(x)) # ガウシアンノイズを生成 plt.figure(figsize=(5,3)) plt.plot(x,data)
フィルタ実行する関数を定義して、上記のdataに適用する。注意点としては、signal.buttordとsignal.butterのアナログオプションをデフォルトのFalse (オプションを指定しない)に戻しておく必要がある。
#バターワースフィルタ def bandpass(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, "band") #フィルタ伝達関数の分子と分母を計算 y = signal.filtfilt(b, a, x) #信号に対してフィルタをかける return y # フィルタ関数を実行 data_bandfilt = bandpass(data, samplerate, fp, fs, gpass, gstop)
2. 信号の周波数解析
次にdataとフィルタを適用したdata_bandfiltの周波数成分を調べる。fft関数はいろいろなモジュールで提供されているが、今回はNumpyのfftモジュールを使用する。
以下のようにnp.fft.fft(data)とするとfftamp0は複素数で数値が返ってくるので、例えば絶対値をとってさらにlogに20をかけることで相対的なdB表示することができる。
（正確にはリファレンスになる電圧やパワーレベルで割っておく必要があるので、ここで使っているdBの絶対値は正しくない）
周波数軸は10秒間のデータなので周波数増分は1/10=0.1Hzとなり、サンプリングレートは2500なので0から0.1Hzずつ2500Hzまでがフルスケールになるが、ナイキスト定理により1250Hzまでが有用な周波数範囲となる。
fftamp0=np.fft.fft(data) fftamp1=np.fft.fft(data_bandfilt) df=0.1 fx=np.arange(0,2500,0.1) np.size(fx) plt.plot(fx,20*np.log10(np.abs(fftamp0))) plt.plot(fx,20*np.log10(np.abs(fftamp1))) plt.xlim([497,503]) plt.vlines(499.5,-40,60,'green','dashed') plt.vlines(500.5,-40,60,'green','dashed') plt.xlim([495,505]) plt.ylim([-40,60]) plt.legend(['Raw','Butterworth band'])
縦の緑破線は通過周波数端の499.5Hzと500.5Hzを示している。バンドパスフィルタ適用後の信号データは通過帯域以外は減衰していることがわかる。
3. バタワース以外の効果
上のバタワースフィルタと同じ手順で第１種チェビシェフ、第２種チェビシェフ、楕円フィルタでバンドパスフィルタを設計して信号データに適用する。 通過、阻止帯域はバタワースと同じだが、阻止域最小損失を各フィルタで特性がひずまないよう程度の値まで大きくしている。以下のコードはバタワースからの主な変更のみを示し、fftやプロットの部分はほぼ胸中なので省略する。
#第１種チェビシェフフィルタ def bandpass_ch1(x, samplerate, fp, fs, gpass, gstop): fn = samplerate / 2 wp = fp / fn ws = fs / fn N, Wn = signal.cheb1ord(wp, ws, gpass, gstop) b, a = signal.cheby1(N,3, Wn, "band") y = signal.filtfilt(b, a, x) return y fp = np.array([499.5,500.5]) # 通過域端周波数[Hz] fs = np.array([499,501]) # 阻止域端周波数[Hz] gpass = 3 # 通過域端最大損失[dB] gstop = 35 # 阻止域端最小損失[dB] data_bandfilt = bandpass_ch1(data, samplerate, fp, fs, gpass, gstop)
第１種チェビシェフフィルタは通過帯域のリップルを許容するため、信号データの周波数成分に対して差が生じている。
#第２種チェビシェフフィルタ def bandpass_ch2(x, samplerate, fp, fs, gpass, gstop): fn = samplerate / 2 wp = fp / fn ws = fs / fn N, Wn = signal.cheb2ord(wp, ws, gpass, gstop) b, a = signal.cheby2(N,gstop, Wn, "band") y = signal.filtfilt(b, a, x) return y fp = np.array([499.5,500.5]) # 通過域端周波数[Hz] fs = np.array([499,501]) # 阻止域端周波数[Hz] gpass = 3 # 通過域端最大損失[dB] gstop = 35 # 阻止域端最小損失[dB] data_bandfilt = bandpass_ch2(data, samplerate, fp, fs, gpass, gstop)
第２種チェビシェフフィルタは通過帯域のリップルはないため、信号データの周波数成分は第１種チェビシェフフィルタよりも一致している。
#楕円フィルタ def bandpass_ell(x, samplerate, fp, fs, gpass, gstop): fn = samplerate / 2 wp = fp / fn ws = fs / fn N, Wn = signal.ellipord(wp, ws, gpass, gstop) b, a = signal.ellip(N,3,gstop, Wn, "band") y = signal.filtfilt(b, a, x) return y fp = np.array([499.5,500.5]) # 通過域端周波数[Hz] fs = np.array([499,501]) # 阻止域端周波数[Hz] gpass = 3 # 通過域端最大損失[dB] gstop = 55 # 阻止域端最小損失[dB] data_bandfilt = bandpass_ell(data, samplerate, fp, fs, gpass, gstop)
楕円フィルタは通過帯域と阻止帯域のリップルを両方許容するかわりに急峻に信号をカットするが、第１種チェビシェフフィルタと同様に通過帯域に対して信号データの周波数成分と差が生じている。。
3. まとめ
今回は信号データの周波数解析とバタワース、チェビシェフ、楕円バンドパスフィルタを適用した場合の周波数成分の変動について紹介した。
いろいろアップデート
今週はiOS15が正式リリースされたが、手持ちのiPad ProもiPad OS15へアップデートしてみた。
まだ、あまり試せていないがiPencilを用いた手書き入力がかなり進化している。
手書きの日本語、漢字をテキストに変換してくれる。ただ数字や記号など崩れていると誤変換が起きやすい。
iPadOSのアップデートに合わせて、手持ちのガジェットのファームウェアもいろいろアップデートしてみた。
まずオーディオテクニカのワイアレスイヤホン
特に不便は感じたことがないが、2年ぐらい前の購入からファームウェアの更新ができるようなのでアプリ(Connect)経由でRev2000からRev2001にアップデート。
片方づつ更新する必要があったりファンクションキーの操作がわかりにくかったり、結構時間がかかったが一応更新に成功した。ファームウェアをアップデートしても特に使い勝手は今まで通り。
次にガーミンのForeathlete 235J。
去年一度ファームウェアをアップデートしたが、今回改めて更新がないかチェックしてみると8.4というバージョンにアップデート可能だった。
こちらはPCにインストールしたGarmin Expressというソフトウェアを用いてUSB接続すれば簡単にアップデートすることができた。こちらも使い勝手は変わらないが、前回更新したときはGPS信号の捕捉が改善されていたので、今回もGPSの改善に期待したい。（前回更新後からだんだんGPS信号捕捉に時間がかかってきている）
