代表的なFIRフィルタである以下フィルタの設計・評価のためのコードをまとめており、フィルタの周波数特性や位相特性、群遅延プロット、信号のフィルタ処理出力ができます。
評価用のサンプル信号について
フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力するコードです。
サイン波にノイズを足したものをサンプル信号とします。
dt = 0.005
fs = 1/dt
times = np.arange(0,10,dt)
N = times.shape[0]
f = 5
sigma = 0.5
np.random.seed(1)
x_s =np.sin(2 * np.pi * times * f)
x = x_s + sigma * np.random.randn(N)
x_np = np.array(x)
y_s = np.zeros(times.shape[0])
y_s[:times.shape[0]//2] = 1
y = y_s + sigma * np.random.randn(N)
plt.figure(figsize=(12,6))
plt.subplot(3,1,1)
plt.plot(times, x, c="red", label=("sin+noise"))
plt.xlim(0,1)
plt.yticks([-2, -1.5, -1, -0.5, 0, 0.5,1, 1.5, 2])
plt.grid(which = "major", axis = "y", color = "gray", alpha = 0.8,
linestyle = "--", linewidth = 1)
plt.legend()
def FFT(Raw_Data, sampling_time):
sampling_time_f = 1/sampling_time
fft_nseg = 256
fft_overlap = fft_nseg // 2
f_raw, t_raw, Sxx_raw = signal.stft(Raw_Data, sampling_time_f, nperseg=fft_nseg, noverlap=fft_overlap)
return f_raw, t_raw, Sxx_raw
def get_start_end_index(t, start_t, end_t):
start_ind=0
while(t[start_ind] < start_t):
start_ind += 1
end_ind = start_ind + 1
while (t[end_ind] < end_t):
end_ind += 1
return start_ind, end_ind
def ave_spectrum(Sxx, t, start_t, end_t):
start_ind, end_ind = get_start_end_index(t, start_t, end_t)
ave_spectrum = np.zeros(Sxx.shape[0])
for i in range(Sxx.shape[0]):
ave_spectrum[i] = np.average(Sxx[i, start_ind:end_ind])
return ave_spectrum
f_raw, t_raw, Sxx_raw = FFT(x_np, dt)
plt.subplot(3,1,2)
plt.pcolormesh(t_raw, f_raw, np.abs(Sxx_raw),vmin=0, vmax=0.4, cmap="jet")
plt.subplot(3,1,3)
start = 0
end = 2
gain = ave_spectrum(np.abs(Sxx_raw), t_raw, start, end)
plt.plot(f_raw, gain, c="blue", label=("normal"))
plt.legend()
plt.show()
前の記事と同じものですが、実行結果は以下の通りです。
グラフ1つ目がサンプル信号の時系列です
グラフ2つ目がサンプル信号の周波数特性を示した時系列です。このグラフは関数「FFT」の出力結果で、横軸:時間 縦軸:周波数 色が出力(db)を示しており、周波数特性を時系列的に見たい場合に有効です。この画像だと5Hz付近で赤色となっており、グラフ3つ目の周波数特性とも合っていることが分かります。
グラフ3つ目がサンプル信号の周波数特性です。「start」「end」に時間を入力し、その時間範囲のデータでの周波数特性を示します。関数「FFT」の結果を関数「get_start_end_index」「ave_spectrum」で切り分け、平均化処理させた結果を描画しています。
Butterworthフィルタについて
代表的なFIRフィルタで、BPF/LPF/HPFと使い勝手が良いです。
通過帯域は全てを通すフィルタであり、リップル(揺らぎ)が無くいため、除去帯域近くの動きも滑らかです。※高次になれば傾斜は急になります
また、通過帯域は今から説明する他のフィルタと比較して線形な位相応答を示しており、周波数ごとの位相遅れの差もほとんど無いのが特徴です。
以下がサンプルコードとなり、フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力する関数です。
「ftype」数字によってLPF/HPF/BPFを選択できるようになっており、設定結果及び「fl」「fh」を任意の周波数/「order」次数を設定することでフィルタのパラメータを設定します。
関数「character_MA」がフィルタ特性の出力用関数
関数「character_Butter」が実際のフィルタ処理用関数になります。
dt = 0.005
fs = 1/dt
ftype = 0
fl = None
fh = 20
order = 4
class Butter_Filter():
def __init__(self, master=None):
print("Butterworse_Filter")
def get_filtertype(self, selected_num):
if selected_num == 0:
ftype1 = "lowpass"
elif selected_num == 1:
ftype1 = "highpass"
else:
ftype1 = "bandpass"
return ftype1
def get_filterparam(self, dt, fl, fh, order, filtertype):
dt = float(dt)
fs = 1 / dt
fn = fs / 2
if filtertype == "bandpass":
Wn = [float(fl)/fn, float(fh)/fn]
N=order//2
elif filtertype == "lowpass":
Wn = float(fh)/fn
N=order
else:
Wn = float(fl)/fn
N=order
return Wn, N
def character_Butter(self, dt, order, fl, fh, ftype):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
b, a = signal.butter(N, Wn, ftype)
w, h = signal.freqz(b, a)
x_freq_rspns = w / np.pi * fn
y_freq_rspns = db(abs(h))
p = np.angle(h) * 180 / np.pi
w, gd = signal.group_delay([b, a])
x_gd = w / np.pi * fn
y_gd = gd
return x_freq_rspns, y_freq_rspns, p, x_gd, y_gd
def Butter_filter(self, data, dt, order, fl, fh, ftype, worN=8192):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
print(Wn, N)
b, a = signal.butter(N, Wn, ftype)
zi = signal.lfilter_zi(b, a)
z, _ = signal.lfilter(b, a, data, zi=zi*data[0])
return z
BFilter = Butter_Filter()
ftype = BFilter.get_filtertype(ftype)
print(ftype)
y = BFilter.Butter_filter(x, dt, order, fl, fh, ftype)
x_freq_rspns_bu, y_freq_rspns_bu, p_bu, x_gd_bu, y_gd_bu = BFilter.character_Butter(dt, order, fl, fh, ftype)
フィルタ特性について
フィルタ特性確認用のサンプルコードです
plt.figure(figsize=(12,6))
plt.subplot(3, 1, 1)
plt.plot(x_freq_rspns_bu, y_freq_rspns_bu)
plt.xlim(0, 100)
plt.ylim(-200, 2)
plt.ylabel('Amplitude [dB]')
plt.grid()
plt.subplot(3, 1, 2)
plt.plot(x_freq_rspns_bu, p_bu)
plt.xlim(0, 100)
plt.ylim(-200, 200)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [deg]')
plt.grid()
plt.subplot(3, 1, 3)
plt.plot(x_gd_bu, y_gd_bu)
plt.xlim(0, 100)
plt.ylim(-5, 15)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Group delay [samples]')
plt.grid()
plt.show()
実行結果は以下の通りです。
※サンプルコードはローパスフィルタ設定
グラフ1つ目がフィルタの周波数特性
グラフ2つ目がフィルタの位相特性
グラフ3つ目がフィルタの群遅延を示しています。
群遅延は、フィルタ処理によって何サンプル分信号が遅れるのかを示しており、その確認ができます。
フィルタ処理結果の確認
フィルタ処理結果確認用のサンプルコードです。
「x」に任意の信号を入れていただくことで処理ができますので、処理したい信号を「x」に入れて試してみてください。
plt.figure(figsize=(12,6))
plt.subplot(3,1,1)
plt.plot(times, x, c="red", label=("sin+noise"))
plt.plot(times, x_s, c="blue", label=("sin"))
plt.plot(times, y, c="green", label=("Butter"))
plt.xlim(0,1)
plt.yticks([-2, -1.5, -1, -0.5, 0, 0.5,1, 1.5, 2])
plt.grid(which = "major", axis = "y", color = "gray", alpha = 0.8,
linestyle = "--", linewidth = 1)
plt.legend()
f_raw, t_raw, Sxx_raw = dfilter.FFT(y, dt)
plt.subplot(3,1,2)
plt.pcolormesh(t_raw, f_raw, np.abs(Sxx_raw),vmin=0, vmax=0.4, cmap="jet")
plt.subplot(3,1,3)
start = 0
end = 2
f_raw, t_raw, Sxx_raw = dfilter.FFT(x, dt)
gain = dfilter.ave_spectrum(np.abs(Sxx_raw), t_raw, start, end)
plt.plot(f_raw, gain, c="red", label=("sin+noise"))
f_raw, t_raw, Sxx_raw = dfilter.FFT(x_s, dt)
gain = dfilter.ave_spectrum(np.abs(Sxx_raw), t_raw, start, end)
plt.plot(f_raw, gain, c="blue", label=("sin"))
f_raw, t_raw, Sxx_raw = dfilter.FFT(y, dt)
gain = dfilter.ave_spectrum(np.abs(Sxx_raw), t_raw, start, end)
plt.plot(f_raw, gain, c="green", label=("Butter"))
plt.legend()
plt.show()
グラフ1つ目の緑プロットがButterworthフィルタ後の出力となります。※フィルタ処理前が赤プロット
グラフ2つ目がフィルタ後出力の周波数特性を示した時系列です。
グラフ3つ目の緑プロットが移動平均フィルタ後の周波数特性になります。
設計通り、20Hz付近から減衰するようなフィルタ特性を持っており、グラフ3つ目の緑プロットと赤プロットを比較して、設計通りのフィルタが作成できていることが分かります。
Chebyshevフィルタについて
通過帯域でのリップルが存在するが、Butterworseと比較して減衰頻度を大きく取れるフィルタ(第一種Chebyshevフィルタ)
減衰頻度をさらに大きく取れるように設計されたフィルタが第二種Chebyshevフィルタ
第一種Chebyshevフィルタ(フィルタ特性/フィルタ処理)
以下がサンプルコードとなり、フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力する関数です。
「ftype」数字によってLPF/HPF/BPFを選択できるようになっており、設定結果及び「fl」「fh」を任意の周波数/「order」次数を設定することでフィルタのパラメータを設定します。
第一種Chebyshevフィルタでは、「attenuation」でリップルについても調整可能です。
dt = 0.005
fs = 1/dt
ftype = 0
fl = None
fh = 20
order = 4
attenuation = 5
class Chebyshev_Filter():
def __init__(self, master=None):
print("Chebyshev_Filter")
def get_filtertype(self, selected_num):
if selected_num == 0:
ftype1 = "lowpass"
elif selected_num == 1:
ftype1 = "highpass"
else:
ftype1 = "bandpass"
return ftype1
def get_filterparam(self, dt, fl, fh, order, filtertype):
dt = float(dt)
fs = 1 / dt
fn = fs / 2
if filtertype == "bandpass":
Wn = [float(fl)/fn, float(fh)/fn]
N=order//2
elif filtertype == "lowpass":
Wn = float(fh)/fn
N=order
else:
Wn = float(fl)/fn
N=order
return Wn, N
def character_Chebyshev(self, dt, order, fl, fh, ripple, ftype):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
b, a = signal.cheby1(N, ripple, Wn, ftype)
w, h = signal.freqz(b, a)
x_freq_rspns = w / np.pi * fn
y_freq_rspns = db(abs(h))
p = np.angle(h) * 180 / np.pi
w, gd = signal.group_delay([b, a])
x_gd = w / np.pi * fn
y_gd = gd
return x_freq_rspns, y_freq_rspns, p, x_gd, y_gd
def Chebyshev_filter(self, data, dt, order, fl, fh, ripple, ftype):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
b, a = signal.cheby1(N, ripple, Wn, ftype)
zi = signal.lfilter_zi(b, a)
z, _ = signal.lfilter(b, a, data, zi=zi*data[0])
return z
CFilter = dfilter.Chebyshev_Filter()
ftype = CFilter.get_filtertype(ftype)
print(ftype)
y = CFilter.Chebyshev_filter(x, dt, order, fl, fh, attenuation, ftype)
x_freq_rspns_ch, y_freq_rspns_ch, p_ch, x_gd_ch, y_gd_ch = CFilter.character_Chebyshev(dt, order, fl, fh, ripple, ftype)
フィルタ特性は以下の通りです。
フィルタ次数にもよりますが、Butterworthフィルタと比較して減衰が大きくなっているように見えます。
ただ、位相/群遅延プロットを見ると、周波数によって遅れが変わってきているので、使う際はこの辺を理解したうえで使う必要ありそうです。
フィルタ処理結果は以下の通りです。
ローパスフィルタとしてうまく機能はしてそうですが、リップルの影響か、sin波の強度が若干低下していますね。
この辺も理解したうえで使う必要ありそうです。
第二種Chebyshevフィルタ(フィルタ特性/フィルタ処理)
第一種Chebyshevフィルタよりも減衰頻度を大きく取れるように設計されたフィルタ
通過帯域にリップルがない代わりに、除去帯域で等リップル性を持つ。
以下がサンプルコードとなり、フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力する関数です。
「ftype」数字によってLPF/HPF/BPFを選択できるようになっており、設定結果及び「fl」「fh」を任意の周波数/「order」次数を設定することでフィルタのパラメータを設定します。
第二種Chebyshevフィルタでは、「attenuation」でリップルについても調整可能です。
dt = 0.005
fs = 1/dt
ftype = 0
fl = None
fh = 20
order = 4
attenuation = 10
class Chebyshev_second_Filter():
def __init__(self, master=None):
print("Chebyshev_Filter")
def get_filtertype(self, selected_num):
if selected_num == 0:
ftype1 = "lowpass"
elif selected_num == 1:
ftype1 = "highpass"
else:
ftype1 = "bandpass"
return ftype1
def get_filterparam(self, dt, fl, fh, order, filtertype):
dt = float(dt)
fs = 1 / dt
fn = fs / 2
if filtertype == "bandpass":
Wn = [float(fl)/fn, float(fh)/fn]
N=order//2
elif filtertype == "lowpass":
Wn = float(fh)/fn
N=order
else:
Wn = float(fl)/fn
N=order
return Wn, N
def character_Chebyshev_second(self, dt, order, fl, fh, attenuation, ftype):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
b, a = signal.cheby2(N, attenuation, Wn, ftype)
w, h = signal.freqz(b, a)
x_freq_rspns = w / np.pi * fn
y_freq_rspns = db(abs(h))
p = np.angle(h) * 180 / np.pi
w, gd = signal.group_delay([b, a])
x_gd = w / np.pi * fn
y_gd = gd
return x_freq_rspns, y_freq_rspns, p, x_gd, y_gd
def Chebyshev_second_filter(self, data, dt, order, fl, fh, attenuation, ftype):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
b, a = signal.cheby2(N, attenuation, Wn, ftype)
zi = signal.lfilter_zi(b, a)
z, _ = signal.lfilter(b, a, data, zi=zi*data[0])
return z
CsFilter = Chebyshev_second_Filter()
ftype = CsFilter.get_filtertype(ftype)
print(ftype)
y = CsFilter.Chebyshev_second_filter(x+1, dt, order, fl, fh, attenuation, ftype)
x_freq_rspns_ch2, y_freq_rspns_ch2, p_ch2, x_gd_ch2, y_gd_ch2 = CsFilter.character_Chebyshev_second(dt, order, fl, fh, attenuation, ftype)
フィルタ特性は以下の通りです。
フィルタ処理結果は以下の通りです。
Besselフィルタについて
最後にBesselフィルタについて紹介します。
Butterworthと非常に似た特性を持つフィルタで、過渡応答性がより滑らかなフィルタ。
通過帯域では線形な位相特性を持つため、信号の遅れも周波数によらずある程度一定を保てることができます。
正直言うとButterworthとBesselフィルタとの使い分けはわかりませんでした。もう少し勉強します。
フィルタ特性/フィルタ処理
dt = 0.005
fs = 1/dt
ftype = 0
fl = None
fh = 20
order = 4
class Bessel_Filter():
def __init__(self, master=None):
print("Chebyshev_Filter")
def get_filtertype(self, selected_num):
if selected_num == 0:
ftype1 = "lowpass"
elif selected_num == 1:
ftype1 = "highpass"
else:
ftype1 = "bandpass"
return ftype1
def get_filterparam(self, dt, fl, fh, order, filtertype):
dt = float(dt)
fs = 1 / dt
fn = fs / 2
if filtertype == "bandpass":
Wn = [float(fl)/fn, float(fh)/fn]
N=order//2
elif filtertype == "lowpass":
Wn = float(fh)/fn
N=order
else:
Wn = float(fl)/fn
N=order
return Wn, N
def character_Bessel(self, dt, order, fl, fh, ftype):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
b, a = signal.bessel(N, Wn, ftype)
w, h = signal.freqz(b, a)
x_freq_rspns = w / np.pi * fn
y_freq_rspns = db(abs(h))
p = np.angle(h) * 180 / np.pi
w, gd = signal.group_delay([b, a])
x_gd = w / np.pi * fn
y_gd = gd
return x_freq_rspns, y_freq_rspns, p, x_gd, y_gd
def Bessel_filter(self, data, dt, order, fl, fh, ftype, worN=8192):
fs = 1 / dt
fn = fs / 2
Wn, N = self.get_filterparam(dt, fl, fh, order, ftype)
b, a = signal.bessel(N, Wn, ftype)
zi = signal.lfilter_zi(b, a)
z, _ = signal.lfilter(b, a, data, zi=zi*data[0])
return z
BeFilter = Bessel_Filter()
ftype = BeFilter.get_filtertype(ftype)
print(ftype)
y = BeFilter.Bessel_filter(x, dt, order, fl, fh, ftype)
x_freq_rspns_be, y_freq_rspns_be, p_be, x_gd_be, y_gd_be = BeFilter.character_Bessel(dt, order, fl, fh, ftype)
フィルタ特性は以下の通りです。
確かにButterworthフィルタと比較して減衰域が滑らかになっているように見えますね。ただ、通過域と減衰域の境が滑らかになっており攻めたフィルタをしたときに境目付近の信号は減衰されてしまう可能性がありそうですね。
フィルタ処理結果は以下の通りです。
