内容
有限インパルス応答Finite Impulse Response (FIR)の設計
無限インパルス応答Infinite Impulse Response (IIR)の設計
はじめに
MNE tutorialのBackground information on filteringを試してみました。
周波数フィルタは時系列データを周波数データに変換して窓関数をかけて特定の周波数を除去し時系列データに戻す方法です。
図示すると下の図のようになり、フーリエ変換・窓関数・逆フーリエ変換部分が周波数フィルタ処理部分になります。
（ホントはフーリエ変換を発展させたラプラス変換の離散版、z変換なのですが、話を簡単にします）
 |
|||
 |
|||
ホームページの縦書きって難しですね。
周波数フィルタは信号雑音比を良くするのに使いますが、肝心な生体データもフィルタの影響で歪んでしまうことがあるので注意しましょう。
周波数フィルタの設計は奥が深いです。本家のホームページと比べてここでは原理的な部分は触れていません。何とか使えそうということが分かれば、とりあえずは十分と思ってます。
Spyderを使います。Python consoleだとdefで関数の定義ができず、IPythonコンソールだとウインドウの制御ができそうもないので、Python consoleを使います。
有限インパルス応答Finite Impulse Response (FIR)の設計
いろいろ導入していきます。
>>>は面倒なので無視します。
import numpy as np
from scipy import signal, fftpack
import matplotlib.pyplot as plt
from mne.time_frequency.tfr import morlet
import mne
時間周波数プロット用の設定です。
sfreq=1000 # サンプリング周波数
f_p=40. # 40Hz
ylim=[-60,10] # 周波数・gain図用 縦軸 dB
xlim=[2,sfreq /2.] # 周波数・gain図用 横軸 Hz
blue='#1f77b4' # たぶん16進数
ローパスフィルタを設計します。
nyq=sfreq/2. # Nyquist周波数=サンプリング周波数÷2
freq=[0,f_p,f_p,nyq] # 窓関数の定義
gain=[1,1,0,0] # 0〜40Hzは1 40Hz以上は0
plt.plot(freq,gain) # こんな感じの窓関数です
横軸は周波数、縦軸はgainです。グラフはこれでいいんですが、関数機能defを使って上図を書くプログラムを書いてみます。
一旦ウインドウは閉じます。
plt.close(1)
図を提示するbox_off関数と図にグラフを書くplot_ideal関数を定義します。
Pythonは行の先頭の場所に意味があります。タブなどで設定する4文字空白は必ず揃えてください。
def box_off(ax):
ax.grid(zorder=0)
for key in ('top','right'):
ax.spines[key].set_visible(False)
:を忘れると改行が・・・とならず、エラーになります。また・・・後に改行して>>>が出るまでが関数定義の部分です。
def plot_ideal(freq,gain,ax):
freq=np.maximum(freq,xlim[0]) # xlim[0]以上の周波数
xs,ys=list(),list() # Matlabだとxs=[];ys=[]
my_freq,my_gain=list(),list()
for ii in range(len(freq)):
xs.append(freq[ii]) # Matlabだとxs=[xs,freq(ii)]
ys.append(ylim[0])
if ii<len(freq)-1 and gain[ii]!=gain[ii+1]:
xs+=[freq[ii],freq[ii+1]] # Matlabだとxs=[xs,freq(ii),freq(ii+1)]
ys+=[ylim[1]]*2
my_freq+=np.linspace(freq[ii],freq[ii+1],20,endpoint=False).tolist()
my_gain+=np.linspace(gain[ii],gain[ii+1],20,endpoint=False).tolist()
else:
my_freq.append(freq[ii])
my_gain.append(gain[ii])
my_gain=10*np.log10(np.maximum(my_gain,10**(ylim[0]/10.))) # dB表示です
ax.fill_between(xs,ylim[0],ys,color='r',alpha=0.1) # 横軸は対数
ax.semilogx(my_freq,my_gain,'r--',alpha=0.5,linewidth=4,zorder=3)
xticks=[1,2,4,10,20,40,100,200,400]
ax.set(xlim=xlim,ylim=ylim,xticks=xticks,xlabel='Frequency (Hz)',ylabel='Amplitude (dB)')
ax.set(xticklabels=xticks)
box_off(ax)
周波数応答を図にしてみます。通常のウインドウの高さを半分にしています。尚、IPythonコンソールだとうまく絵が描けません。
half_height=np.array(plt.rcParams['figure.figsize'])*[1,0.5] # ウィンドウの高さを半分に
ax=plt.subplots(1,figsize=half_height)[1]
plot_ideal(freq,gain,ax) # 先ほど定義した関数です
ax.set(title='Ideal %s Hz lowpass' % f_p)
mne.viz.tight_layout()
plt.show()
ウインドウを閉じてないとFigure 2となります。これが目的とするフィルタです。矩形の窓関数となっています。これをフーリエ変換すると矩形波はsinc関数になります。sinc関数は−∞〜＋∞の大きさがあります。デジタルですんで無限大の大きさのsinc関数を設定するにはいきませんので、以下のようにします。
n=int(round(0.1*sfreq))+1 # sinc関数の大きさ
t=np.arange(-n//2,n//2)/sfreq # sinc関数をπに変換
h=np.sinc(2*f_p*t)/(4*np.pi) # 作成したフィルタとして用いるsinc関数
n
t
h
>>>と反応が時々ずれることがありますが、これは愛嬌です。
sinc関数の大きさnは101です。このフィルタの周波数・gain図（周波数応答）を見てみます。
plot_filter関数を定義します。
def plot_filter(h,title,freq,gain,show=True):
if h.ndim==2: # 2次セクション(項)型の場合
sos=h # Second Order Sectionの頭文字
n=mne.filter.estimate_ringing_samples(sos)
h=np.zeros(n)
h[0]=1
h=signal.sosfilt(sos,h) # 2次セクション型を伝達関数型に変換
H=np.ones(512,np.complex128)
for section in sos:
f,this_H=signal.freqz(section[:3],section[3:])
H*=this_H # 2次セクション
else: # 伝達関数型の場合
f,H=signal.freqz(h)
fig,axs=plt.subplots(2)
t=np.arange(len(h))/sfreq
axs[0].plot(t,h,color=blue)
axs[0].set(xlim=t[[0,-1]],xlabel='Time (sec)',ylabel='Amplitude h(n)',title=title)
box_off(axs[0])
f*=sfreq/(2*np.pi)
axs[1].semilogx(f,10*np.log10((H*H.conj()).real),color=blue,linewidth=2,zorder=4)
plot_ideal(freq,gain,axs[1])
mne.viz.tight_layout()
if show:
plt.show()
実際に図示してみます。
plot_filter(h,'Sinc (0.1 sec)',freq,gain)
上手はsinc関数、下図が窓関数＝周波数・gain図＝周波数応答です。あまりシャープなフィルタにはなってません。sinc関数の幅を広げ、次数を約10倍にしてフィルタを設計してみます。
n=int(round(1.*sfreq))+1
t=np.arange(-n//2,n//2)/sfreq
h=np.sinc(2*f_p*t)/(4*np.pi)
plot_filter(h,'Sinc(1.0 sec)',freq,gain)
sinc関数の幅は10倍に広がってます。周波数応答は前よりいい感じになってます。
次数を約100倍にしてみます。
n=int(round(10.*sfreq))+1
t=np.arange(-n//2,n//2)/sfreq
h=np.sinc(2*f_p*t)/(4*np.pi)
plot_filter(h,'Sinc(10.0 sec)',freq,gain)
sinc関数の幅は100倍になっています。周波数応答はさらにいい感じになってます。このフィルタ、データサンプルが10001ないと使えません。
sinc関数は周波数フィルタの基本中の基本なんですが、このまま使うことはありません。MATLABにはいろいろな関数が用意されていて、そのうち以下のものはscipyでも用意されています。
1. scipy.signa.remz() MATLABのfirpm関数 Parks-McClellan最適等リップルFIRフィルタ
2. scipy.signal.firwin2() MATLABのfir2関数 周波数サンプリング手法を使ったFIRの任意整形フィルタ
3. scipy.signal.firls() MATLABのfirls関数 最小二乗線形位相FIRフィルタ
4. その他 Frequency domain design (construct filter in Fourier domain and use an IFFT to invert it)
有限の大きさのフィルタで厳密な周波数応答を求めるのは諦めます。移行帯を設けた窓関数でフィルタの設計をすることにします。以下のような周波数応答を目標とします。
trans_bandwidth=10 # 10㎐移行帯を設ける
f_s=f_p+trans_bandwidth # 40~50Hz
freq=[0.,f_p,f_s,nyq]
gain=[1.,1.,0.,0.]
ax=plt.subplots(1,figsize=half_height)[1]
plot_ideal(freq,gain,ax)
ax.set(title='%s Hz lowpass width a %s Hz transition' % (f_p, trans_bandwidth))
mne.viz.tight_layout()
plt.show()
上述のscipy.signal.firwin2()を使ってフィルタを設計します。
h=signal.firwin2(n,freq,gain,nyq=nyq)
plot_filter(h,'Windowed 10-Hz transition (10.0 sec)',freq,gain)
フィルタの次数が高いのでこれくらいは当たり前です。次数を減らします。
n=int(round(sfreq*0.5))+1
h=signal.firwin2(n,freq,gain,nyq=nyq)
plot_filter(h,'Windowed 10-Hz transition (0.5 sec)',freq,gain)
もっと次数を減らします。
n=int(round(sfreq*0.2))+1
h=signal.firwin2(n,freq,gain,nyq=nyq)
plot_filter(h,'Windowed 10-Hz transition (0.5 sec)',freq,gain)
移行帯を25Hzにしてみます。
trans_bandwidth=25
f_s=f_p+trans_bandwidth
freq=[0,f_p,f_s,nyq]
h=signal.firwin2(n,freq,gain,nyq=nyq)
plot_filter(h,'Windowed 50-Hz transition (0.2 sec)',freq,gain)
煩雑なのでウインドウは全部閉じます。
for ii in [1,2,3,4,5,6,7,8,9]:
plt.close(ii) # ウインドウを閉じる
有限インパルス応答の適応
ガウス窓の正弦曲線（例Morlet waveletの虚部）に白色雑音などを足したデータで検討します。正弦曲線部分の周波数帯域はローパスフィルタの周波数帯域に含まれるとします。
Morletのwaveletの準備をします。
dur=10. # 10秒
center=2.
morlet_freq=f_p # 40.0
tlim=[center-0.2,center+0.2]
tticks=[tlim[0],center,tlim[1]]
flim=[20,70] # 20〜70Hz
Morletのwaveletを作成します。
x=np.zeros(int(sfreq*dur)) # 1000x10
blip=morlet(sfreq,[morlet_freq],n_cycles=7)[0].imag/20.
n_onset=int(center*sfreq)-len(blip)//2
x[n_onset:n_onset+len(blip)]+=blip
x_orig=x.copy()
雑音を加えます。
rng=np.random.RandomState(0)
x+=rng.randn(len(x))/1000.
x+=np.sin(2.*np.pi*60*np.arange(len(x))/sfreq)/2000.
scipyのsignal.firwin2を使ってFIRフィルタを設計します。
transition_band=0.25*f_p
f_s=f_p+transition_band
filter_dur=6.6/transition_band # 秒
n=int(sfreq*filter_dur)
freq=[0.,f_p,f_s,sfreq/2.]
gain=[1.,1.,0.,0.]
h=signal.firwin2(n,freq,gain,nyq=sfreq/2.)
x_shallow=np.convolve(h,x)[len(h)//2:]
窓関数を見てみます。
plot_filter(h,'MNE-Python 0.14 default',freq,gain)
このフィルタはMNE-Python 0.14の初期設定で使われているものです。
MNE-Python 0.13の初期状態では以下のフィルタを使用していました。
transition_band=0.5 # Hz
f_s=f_p+transition_band
filter_dur=10.0 # 秒
n=int(sfreq*filter_dur)
freq=[0.,f_p,f_s,sfreq/2.]
gain=[1.,1.,0.,0.]
h=signal.firwin2(n,freq,gain,nyq=sfreq/2.)
x_steep=np.convolve(np.convolve(h,x)[::-1],h)[::-1][len(h)-1:-len(h)-1]
plot_filter(h,'MNE-Python 0.13 default',freq,gain)
次数が大きければこうなりますようなぁ。
次はMNE-Cのフィルタです。
h=mne.filter.design_mne_c_filter(sfreq,l_freq=None,h_freq=f_p+2.5)
x_mne_c=np.convolve(h,x)[len(h)//2:]
transition_band=5 # Hz
f_s=f_p+transition_band
freq=[0.,f_p,f_s,sfreq/2.]
gain=[1.,1.,0.,0.]
plot_filter(h,'NME-C default',freq,gain)
元のMNE tutorialのホームページの意訳です。
MNE-Python 0.13とMNE-Cのローパスフィルタはいい感じで高周波をカットしていますが、移行帯でringing（さざ波雑音が長時間続く）が起こりやすい可能性があります。MNE-Pytthon 0.14はringing予防のためMNE-Python 0.14やMNE-Cよりもローパスフィルタはあえて緩くものを使ってます。
まずウインドウを出して関数plot_signalの定義をします。
axs=plt.subplots(1,2)[1]
def plot_signal(x,offset):
t=np.arange(len(x))/sfreq
axs[0].plot(t,x+offset)
axs[0].set(xlabel='Time (sec)',xlim=t[[0,-1]])
box_off(axs[0])
X=fftpack.fft(x)
freqs=fftpack.fftfreq(len(x),1./sfreq)
mask=freqs>=0
X=X[mask]
freqs=freqs[mask]
axs[1].plot(freqs,20*np.log10(np.abs(X)))
axs[1].set(xlim=xlim)
yticks=np.arange(5)/-30.
yticklabels=['Original','Noisy','FIR-shallow (0.14)','FIR-steep (0.13)','FIR-steep (MNE-C)']
で、これを使ってフィルタ処理の結果を示します。
MNE-Python 0.13とMNE-cではさざ波が出ています。
ウインドウは全部閉じることにします。
for ii in [1,2,3,4]:
plt.close(ii) # ウインドウを閉じる
無限インパルス応答Infinite Impulse Response (IIR)の設計
MNE-Pythonはscipy.signalのscipy.signal.iirfilter()やscipy.signal.iirdesign()を使って無限インパルス応答のフィルタを設計しています。Buterworthフィルタを作成してみます。
sos=signal.iirfilter(2,f_p/nyq,btype='low',ftype='butter',output='sos')
plot_filter(sos,'Butterworth order=2',freq,gain)
但し実際には位相のずれを戻すためsosfilffilt()を使っています。
sosfiltfilt=mne.fixes.get_sosfiltfilt()
x_shallow=sosfiltfilt(sos,x)
Butterworthの次数を増やしてみます。
sos=signal.iirfilter(8,f_p/nyq,btype='low',ftype='butter',output='sos')
plot_filter(sos,'Butterworth order=8',freq,gain)
x_shteep=sosfiltfilt(sos,x)
次はChebychevのtype 1です。周波数上のrippleは1dB以下とします。
sos=signal.iirfilter(8,f_p/nyq,btype='low',ftype='cheby1',output='sos',rp=1) # 最大ripple dB
plot_filter(sos,'Chebychev-1 order=8, ripple=1 dB',freq,gain)
x_ripple1=sosfiltfilt(sos,x)
次はChebychevのtype 1で、周波数上のrippleは6dB以下とします。
sos=signal.iirfilter(8,f_p/nyq,btype='low',ftype='cheby1',output='sos',rp=1) # 最大ripple dB
plot_filter(sos,'Chebychev-1 order=8, ripple=6 dB',freq,gain)
x_ripple6=sosfiltfilt(sos,x)
無限インパルス応答の適応
有限インパルス応答で作成したMorletのwavelet信号＋雑音データを使います。
axs=plt.subplots(1,2)[1]
ytcks=np.arange(4)/-30.
yticklabels=['Original','Noisy','Butterworth-2','Butterworth-8']
plot_signal(x_orig,offset=yticks[0])
plot_signal(x,offset=yticks[1])
plot_signal(x_shallow,offset=yticks[2])
plot_signal(x_steep,offset=yticks[3])
axs[0].set(xlim=tlim,title='IIR, Lowpass=%d Hz' % f_p,xticks=tticks,ylim=[-0.125,0.025],yticks=yticks,yticklabels=yticklabels)
for text in axs[0].get_yticklabels():
text.set(rotation=45,size=8)
axs[1].set(xlim=flim,ylim=ylim,xlabel='Frequency (Hz)',ylabel='Magnitude (dB)')
box_off(axs[0])
box_off(axs[1])
mne.viz.tight_layout()
plt.show()
Chebychev1の結果は以下の通りです。
axs=plt.subplots(1,2)[1]
ytcks=np.arange(4)/-30.
yticklabels=['Original','Noisy','Chebychev1-1','Chebychev1-6']
plot_signal(x_orig,offset=yticks[0])
plot_signal(x,offset=yticks[1])
plot_signal(x_ripple1,offset=yticks[2])
plot_signal(x_ripple6,offset=yticks[3])
axs[0].set(xlim=tlim,title='IIR, Lowpass=%d Hz' % f_p,xticks=tticks,ylim=[-0.125,0.025],yticks=yticks,yticklabels=yticklabels)
for text in axs[0].get_yticklabels():
text.set(rotation=45,size=8)
axs[1].set(xlim=flim,ylim=ylim,xlabel='Frequency (Hz)',ylabel='Magnitude (dB)')
box_off(axs[0])
box_off(axs[1])
mne.viz.tight_layout()
plt.show()
このフィルタ、雑音以外の信号も無くなり、かつフィルタの雑音が乗ってしまってます。ダメですね。
ウインドウは全部閉じます。
for ii in [1,2,3,4,5,6]:
plt.close(ii) # ウインドウを閉じる
フィルタ時の留意点
本家のホームページに、フィルタ処理後に新たな信号が生じてしまうとか、ringingやrippleなどいろいろ書いてあります。
元の信号とフィルタ処理後の信号を見比べる習慣は必要だと思います。
ローパスフィルタとハイパスフィルタの影響例
まず信号を作成します。
x=np.zeros(int(2*sfreq))
t=np.arange(0,len(x))/sfreq-0.2
onset=np.where(t>=0.5)[0][0]
cos_t=np.arange(0,int(sfreq*0.8))/sfreq
sig=2.5-2.5*np.cos(2*np.pi*(1./0.8)*cos_t)
x[onset:onset+len(sig)]=sig
次にIIRフィルタの設計です。
iir_lp_30=signal.iirfilter(2,30./sfreq,btype='lowpass')
iir_hp_p1=signal.iirfilter(2,0.1/sfreq,btype='highpass')
iir_lp_2 =signal.iirfilter(2,2. /sfreq,btype='lowpass')
iir_hp_2 =signal.iirfilter(2,2. /sfreq,btype='highpass')
IIRフィルタ処理後、位相のズレがゼロになるように補正したデータを計算します。
x_lp_30=signal.filtfilt(iir_lp_30[0],iir_lp_30[1],x,padlen=0)
x_hp_p1=signal.filtfilt(iir_hp_p1[0],iir_hp_p1[1],x,padlen=0)
x_lp_2 =signal.filtfilt(iir_lp_2[0], iir_lp_2[1], x,padlen=0)
x_hp_2 =signal.filtfilt(iir_hp_2[0], iir_hp_2[1], x,padlen=0)
結果を示します。
xlim=t[[0,-1]] # Matlabだとt(1,end)’
ylim=[-2,6]
xlabel='Time (sec)'
ylabel='Amplitude ($\mu$V)' # Latex!
tticks=[0.,0.5,1.3,t[-1]]
axs=plt.subplots(2,2)[1].ravel()
for ax,x_f,title in zip(axs,[x_lp_2,x_lp_30,x_hp_2,x_hp_p1],['LP$_2$','LP$_{30}$','HP$_2$','LP$_{0.1}$']):
ax.plot(t,x,color='0.5')
ax.plot(t,x_f,color='k',linestyle='--')
ax.set(ylim=ylim,xlim=xlim,xticks=tticks,title=title,xlabel=xlabel,ylabel=ylabel)
box_off(ax)
mne.viz.tight_layout()
plt.show()
実線の元の信号がフィルタ処理後、破線で表示されています。
基線補正の有用性の例
ハイパスフィルターであっても基線の補正はした方が、信号の歪みは少ないとされます。基線の補正というのはある時間帯の平均値分を引き算するということです。
関数を定義して図示します。
def baseline_plot(x):
all_axs=plt.subplots(3,2)[1]
for ri,(axs,freq) in enumerate(zip(all_axs,[0.1,0.3,0.5])):
for ci,ax in enumerate(axs):
if ci==0:
iir_hp=signal.iirfilter(4,freq/sfreq,btype='highpass',output='sos')
x_hp=sosfiltfilt(iir_hp,x,padlen=0)
else:
x_hp-=x_hp[t<0].mean() # ここで基線の補正
ax.plot(t,x,color='0.5')
ax.plot(t,x_hp,color='k',linestyle='--')
if ri==0:
ax.set(title=('No ' if ci==0 else '')+'Baseline Correction')
box_off(ax)
ax.set(xticks=tticks,ylim=ylim,xlim=xlim,xlabel=xlabel)
ax.set_ylabel('%0.1f Hz' % freq,rotation=0,horizontalalignment='right')
mne.viz.tight_layout()
plt.suptitle(title)
baseline_plot(x)
タイトルのLP0.1はFigure1のものがそのまま残ったものです。左が基線補正無し、右がゼロ秒以前のデータで基線を補正したデータです。0.1Hz〜、0.3Hz〜、0.5Hz〜のハイパスフィルター処理前後を表示しています。
ゼロ秒以前の基線補正時間帯に活動があった場合どうなるか、について検討します。
n_pre=(t<0).sum()
sig_pre=1-np.cos(2*np.pi*np.arange(n_pre)/(0.5*n_pre))
x[:n_pre]+=sig_pre
baseline_plot(x)
基線補正で使う時間帯に信号があるので、その平均値分基線が下がってます。また基線が下がった分、フィルタ処理後の波形も変化してます。要するに基線補正で設定する時間帯は定常状態、変な信号を含まないこと！ということですね。
各アプリケーションのフィルタ
MNE-Pythonのフィルタ
mne.io.Raw.filter()という関数を用意してますが、これはmne.filter.filter_data()という関数を使っていて、さらにmne.filter.filter_data()はscipy.signal.firwin2()と位相のズレ補正（ゼロ位相）処理を組み合わせたものを使っています。
変数としてローパスフィルターはl_freq、ハイパスフィルターはh_freq、移行帯はl_trans_bandwidth、h_trans_bandwidthを使うこととしています。l_trans_bandwidth=’auto’、h_trans_bandwidth=’auto’としたら以下の値が適応されます。
|
l_freq h_freqの値 |
l_trans_bandwidth |
h_trans_bandwidth |
|
0.01 |
0.01 |
2.0 |
|
0.1 |
0.1 |
2.0 |
|
1.0 |
1.0 |
2.0 |
|
2.0 |
2.0 |
2.0 |
|
4.0 |
2.0 |
2.0 |
|
8.0 |
2.0 |
2.0 |
|
10.0 |
2.5 |
2.5 |
|
20.0 |
5.0 |
5.0 |
|
40.0 |
10.0 |
10.0 |
|
45.0 |
11.25 |
5.0 |
|
48.0 |
12.0 |
2.0 |
filter_length=’auto’の時、フィルタの長さは最小移行帯の逆数に設定されます。引数としてfir_windowHでann窓、Hamming窓、Blackman窓が指定されば場合、フィルタ長はこの逆数に6.2、6.6、11.0をかけた長さとなります。
MNE-Cのフィルタ
l ローパスフィルタの移行帯は5Hzです
l ハイパスフィルタの移行帯は3サンプル/サンプリング周波数 Hzです。
l フィルタ長は8197サンプルです
本家のホームページには他にもいろいろ書いてあります。
EEGLAB
MNE-Python 0.14の初期設定のフィルタと同じようなフィルタです。
Fieldtrip
初期設定では周波数通過フィルタは4次のButterworth、ローパスは2次のButterworthで、ゼロ位相処理したものです。
MNE-Pythonだとraw.filter(…,method=’iir’)です。詳しくはmne.filter.construct_iir_filter()のoptionの説明を見てください。