バンドパスフィルタについて

懐かしいエントリへのアクセスがあったので.

バンドパスフィルタについて – 教えてgoo!

ちゃんとした話は,ベストアンサーの方や質問者さんの補足に譲るとして.

当該エントリは元々別のブログで書いたもので,移転時にフォーマットが崩れていて見難いし,確かMatlabやらScilabやらを使っていた様な気がするけど,今はそういう環境が手元に無いので,リンクは無し(教えてgooのリンクも無効っぽいので,今更過去のエントリは掘り返さない).

簡単に云えば,移動平均の周波数特性についてのエントリで,

積分⇔加算(移動平均フィルタ;(x[k]+x[k+1])/2)⇔ローパス
微分⇔減算(移動差分フィルタ;(x[k]-x[k+1])/2)⇔ハイパス

一つ離れた隣り合う値の差分を取る→バンドパス・フィルタ(Band Pass Filter;BPF)
一つ離れた隣り合う値の和を取る→バンド・エリミネート・フィルタ(BEF;Band Eliminate Filter,ノッチ・フィルタ;Notch Filter)

矢印の形に注意

と云う話だった.これは,シンプルに考えれば物凄く分かり易い話で(実は別のエントリで補足説明を書いていたりしたんだけど,),

差分を取る→微分(変化を強調)→高周波域が優位→ハイパス
和を取る→積分(平均化され変化が消される)→低周波域が優位→ローパス
で,さらに1つ飛ばしで差分・和を取った場合,
その間の値がそれぞれ強調されたり,平均化されたり.

その結果として,間だけが出力される(BPF),間だけカットされる(BEF)ような利得特性を持つ.
積分回路とハイパス・フィルタ,微分回路とローパス・フィルタが等価ではないように,
これらもその周波数応答,特性が似ているのであって,同等ではないけども,
しかし,これらを理解するうえで,なかなか大事な話かな~と思う.的な話だった.

ちょっと周波数領域といった小難しい概念は捨てて,時間と空間にだけ目を向けてあげると良い.時間~空間なので,時間的拡がりというちょっと難しい世界を,空間的広がりに置き換えてやる.そして,データの中心をゼロに合わせて,ゼロ点から空間的広がりを持つデータという風に考えると,バンドパスというのは「任意の範囲のデータが欲しい」という事で,バンドエリミネートは「任意の範囲のデータを消したい」という事なので,ゼロを中心にある範囲のデータが欲しい場合は,ゼロを中心に合わせて[…, 0, 0, 1, 1, 1, 0, 0, …]みたいにフィルタリングしてやればいいし,ゼロを中心にある範囲のデータを消したい場合は,ゼロを中心に合わせて[…, 1, 1, 0, 0, 0, 1, 1, …]ってウェイトを掛けてやれば良い.そして,これをみたら気付く人は気付くと思うけど,この考え方はまさにFFTフィルタのそれで,空間と周波数領域の違いはあるけど,考え方は同じ.そして,もう少し実際的な話に考えを戻せば,拡がりを減じる方向の処理というのは「差分」で,拡がりを強調する方向の処理というのは「和」って事だ.

幾つかの実例をみていくと(簡単に済ますけど),

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal


t = np.linspace(0, 1, 10000)#浮動小数点誤差の影響を除くため大きめの値
y = np.sin(2.0*np.pi*t * 10) + np.sin(2.0*np.pi*t * 50) + np.sin(2.0*np.pi*t * 100)

fig = plt.figure(figsize=(13, 10))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

ax1.plot(t, y)
ax1.set_xlabel('Time [sec]')
ax1.set_ylabel('Amplitude')

freq1, P1 = signal.periodogram(y, 10000)
ax2.plot(freq1[:200], 10*np.log10(P1[:200]), 'b', label='periodogram')
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Power/frequency [dB/Hz]')

FireShot Capture 218 - JupyterLab Alpha Preview - http___localhost_8888_lab

10, 50, 100 [Hz]のsin混合波を考える.

[40, 60] [Hz]のバンドパスフィルタを掛けると,

b, a = signal.butter(3, [40/5000, 60/5000], 'bandpass')
fil_y = signal.filtfilt(b, a, y)

fig = plt.figure(figsize=(13, 10))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

ax1.plot(t, fil_y)
ax1.set_xlabel('Time [sec]')
ax1.set_ylabel('Amplitude')

freq1, P1 = signal.periodogram(fil_y, 10000)
ax2.plot(freq1[:200], 10*np.log10(P1[:200]), 'b', label='periodogram')
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Power/frequency [dB/Hz]')

FireShot Capture 219 - JupyterLab Alpha Preview - http___localhost_8888_lab

当然,50 [Hz]だけが残る.この時,フィルタの次数はというと,

b, a =

(array([ 2.44962276e-07, 0.00000000e+00, -7.34886829e-07,
0.00000000e+00, 7.34886829e-07, 0.00000000e+00,
-2.44962276e-07]),
array([ 1. , -5.97203707, 14.86338053, -19.73309648,
14.73937996, -5.87280724, 0.9751803 ]))

FIRフィルタでもやってみると,

a6 = 1
numtaps = 1023
b6 = signal.firwin(numtaps, [40/5000, 60/5000], pass_zero=False)
fir_y = signal.lfilter(b6, a6, y)
delay = (numtaps-1) / 2 * (1/10000)

fig = plt.figure(figsize=(13, 10))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

ax1.plot(t-delay, fir_y)
ax1.set_xlabel('Time [sec]')
ax1.set_ylabel('Amplitude')

freq1, P1 = signal.periodogram(fir_y, 10000)
ax2.plot(freq1[:200], 10*np.log10(P1[:200]), 'b', label='periodogram')
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Power/frequency [dB/Hz]')

FireShot Capture 220 - JupyterLab Alpha Preview - http___localhost_8888_lab

フィルタは

b, a =

(array([ 0.03567007, 0.24113472, 0.44676917, 0.24113472, 0.03567007]), 1)

勿論,これらはバンドパスフィルタだけど,「一つ離れた隣り合う値の差分を取る」形には必ずしもなっていない.あくまで,必要条件ではない(逆命題は必ずしも真では無い).

幾つかみていくと,

b7 = signal.firwin2(5, [0.0, 40/5000, 60/5000, 1.0], [0, 1.0, 1.0, 0.0])
a7 = np.ones(b7.size)

w, h = signal.freqs(b7, a7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency')
plt.ylabel('Amplitude response [dB]')
plt.grid()
print(b7)

[-0.00506073 0.0780304 0.44281377 0.0780304 -0.00506073]

以下,5,7,9・・・

[-0.00301124 -0.01961032 0.11126557 0.44281377 0.11126557 -0.01961032
-0.00301124]

[-0.00253036 -0.00175429 -0.01707996 0.15067301 0.47444332 0.15067301
-0.01707996 -0.00175429 -0.00253036]

[-0.00181868 -0.00530909 -0.00325034 -0.02157603 0.1588362 0.47444332
0.1588362 -0.02157603 -0.00325034 -0.00530909 -0.00181868]

[-0.00253036 -0.0032197 -0.00980516 -0.00441165 -0.02435476 0.16340268
0.47444332 0.16340268 -0.02435476 -0.00441165 -0.00980516 -0.0032197
-0.00253036]

b7 = np.array([-1, 0, 1, 0, -1])
a7 = np.ones(b7.size)

w, h = signal.freqs(b7, a7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency')
plt.ylabel('Amplitude response [dB]')
plt.grid()
print(b7)

[-1 0 1 0 -1]

FireShot Capture 222 - JupyterLab Alpha Preview - http___localhost_8888_lab

b7 = np.array([1, 0, 1])
a7 = np.ones(b7.size)

w, h = signal.freqs(b7, a7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency')
plt.ylabel('Amplitude response [dB]')
plt.grid()
print(b7)

FireShot Capture 223 - JupyterLab Alpha Preview - http___localhost_8888_lab

b7 = np.array([1, 0])
a7 = np.ones(b7.size) * b7.size

w, h = signal.freqs(b7, a7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency')
plt.ylabel('Amplitude response [dB]')
plt.grid()
print(b7)

FireShot Capture 224 - JupyterLab Alpha Preview - http___localhost_8888_lab

b7 = np.array([0, 1])
a7 = np.ones(b7.size) * b7.size

w, h = signal.freqs(b7, a7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency')
plt.ylabel('Amplitude response [dB]')
plt.grid()
print(b7)

FireShot Capture 225 - JupyterLab Alpha Preview - http___localhost_8888_lab
 
 
 
 
 
その他(グラフ):

移動平均はローパス.n数が増えればカットオフ周波数が低域に遷移.

import scipy.fftpack as fft


def _freqz(b, a=[1]):
    N = 5001
    w = np.arange(5001)
    h = fft.fft(b, 2*N)[:N] / fft.fft(a, 2*N)[:N]
    return w, h


b7 = np.array([1, 1])/2

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()
print(b7)

FireShot Capture 230 - JupyterLab Alpha Preview - http___localhost_8888_lab.png

b7 = np.array([1, 1, 1])/3

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()
print(b7)

FireShot Capture 227 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 5
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 231 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 21
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 232 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 25
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 233 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 51
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 234 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 75
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 235 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 121
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 236 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 200
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 237 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 511
b7 = np.ones(N)/N

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 238 - JupyterLab Alpha Preview - http___localhost_8888_lab
 
 
 

・移動差分はハイパス.ただし,差分点を無限小方向に増やさないといけない(マクロにみていけば,結局はマクロな傾向(低周期)しかみえなくなっていく).

N = 2
b7 = np.concatenate([np.array([1]), -np.ones(N-1)])/(N)

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 239 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 3
b7 = np.concatenate([np.array([1]), -np.ones(N-1)])/(N)

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 240 - JupyterLab Alpha Preview - http___localhost_8888_lab

N = 5
b7 = np.concatenate([np.array([1]), -np.ones(N-1)])/(N)

w, h = _freqz(b7)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 241 - JupyterLab Alpha Preview - http___localhost_8888_lab

移動平均は割と色々な場面で用いられる理由,移動差分なんて余りお目にかけない理由がよく分かる.
また,金融市場などでよく使われるSMA(n=25, 120,200)等は,前後(n=21とか51,75,100)と比べても,ちょっとしたn数差でありながら割と綺麗な周波数特性で,多くの人が経験的に分かっていて,その辺りが好まれているのだろうなと推察される.
 
 
 

2つを合体.

b7 = np.array([1, 1])/2
b7_2 = np.array([1, -1])/2


w, h = _freqz(b7)
w2, h2 = _freqz(b7_2)
plt.semilogx(w, 20*np.log10(abs(h*h2)))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.grid()

FireShot Capture 242 - JupyterLab Alpha Preview - http___localhost_8888_lab

広告
カテゴリー: 未分類 パーマリンク

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中