つれづれなる備忘録

日々の発見をあるがままに綴る

Pythonによるデータ処理8 ~ フィルタ処理5

 今回は前回紹介したバタワース、チェビシェフ、楕円特性フィルタを信号データに適用した際の周波数数成分の変化について紹介する。

atatat.hatenablog.com

1. 準備

設計したフィルタを信号データに適用する方法は以前紹介したが、今回はFFTを用いてフィルタ適用前後の周波数成分がどのように変化するか調べる。

atatat.hatenablog.com

まず、バタワースでバンドパスフィルタを設計して信号データに適用する

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.buttordsignal.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. バタワース以外の効果

 上のバタワースフィルタと同じ手順で第1種チェビシェフ、第2種チェビシェフ、楕円フィルタでバンドパスフィルタを設計して信号データに適用する。 通過、阻止帯域はバタワースと同じだが、阻止域最小損失を各フィルタで特性がひずまないよう程度の値まで大きくしている。以下のコードはバタワースからの主な変更のみを示し、fftやプロットの部分はほぼ胸中なので省略する。

#第1種チェビシェフフィルタ
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)

"第1種チェビシェフバンドパスフィルタ"
第1種チェビシェフバンドパスフィルタ

第1種チェビシェフフィルタは通過帯域のリップルを許容するため、信号データの周波数成分に対して差が生じている。

"第1種チェビシェフバンドパスの周波数変動"
第1種チェビシェフバンドパスの周波数変動

#第2種チェビシェフフィルタ
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)

"第2種チェビシェフバンドパスフィルタ"
第2種チェビシェフバンドパスフィルタ

第2種チェビシェフフィルタは通過帯域のリップルはないため、信号データの周波数成分は第1種チェビシェフフィルタよりも一致している。

"第2種チェビシェフバンドパスの周波数変動"
第2種チェビシェフバンドパスの周波数変動

#楕円フィルタ
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)

"楕円バンドパスフィルタ"
楕円バンドパスフィルタ

楕円フィルタは通過帯域と阻止帯域のリップルを両方許容するかわりに急峻に信号をカットするが、第1種チェビシェフフィルタと同様に通過帯域に対して信号データの周波数成分と差が生じている。。

"楕円バンドパスの周波数変動"
楕円バンドパスの周波数変動

3. まとめ

 今回は信号データの周波数解析とバタワース、チェビシェフ、楕円バンドパスフィルタを適用した場合の周波数成分の変動について紹介した。