まるっとワーク

データ分析・電子工作など気になることを残していきます

Python:ディジタルフィルタの設計② Butterworth/Chebyshev/Bessel

Pythonから利用できる数学処理ライブラリ、Scipyを使ったディジタルフィルタの設計・評価方法の続きをまとめていきます。
前回同様、数学的な話やディジタルフィルタについての概要は飛ばし、フィルタ設計・評価に使えるサンプルプログラムを載せただけですが、何かに使えれば幸いです。

前の記事:  
dango-study.hatenablog.jp

本記事で紹介したコード:
github.com


サンプルプログラム

代表的なIIRフィルタである以下フィルタの設計・評価のためのコードをまとめており、フィルタの周波数特性や位相特性、群遅延プロット、信号のフィルタ処理出力ができます。

評価用のサンプル信号について

フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力するコードです。
サイン波にノイズを足したものをサンプル信号とします。

dt = 0.005 #1stepの時間[sec]
fs = 1/dt
times  =  np.arange(0,10,dt)
N = times.shape[0]

f  = 5  #サイン波の周波数[Hz]
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()




前の記事と同じものですが、実行結果は以下の通りです。
f:id:toku_dango:20210612231052p:plain
グラフ1つ目がサンプル信号の時系列です
グラフ2つ目がサンプル信号の周波数特性を示した時系列です。このグラフは関数「FFT」の出力結果で、横軸:時間 縦軸:周波数 色が出力(db)を示しており、周波数特性を時系列的に見たい場合に有効です。この画像だと5Hz付近で赤色となっており、グラフ3つ目の周波数特性とも合っていることが分かります。
グラフ3つ目がサンプル信号の周波数特性です。「start」「end」に時間を入力し、その時間範囲のデータでの周波数特性を示します。関数「FFT」の結果を関数「get_start_end_index」「ave_spectrum」で切り分け、平均化処理させた結果を描画しています。

Butterworthフィルタについて

代表的なフィルタで、BPF/LPF/HPFと使い勝手が良いです。
通過帯域は全てを通すフィルタであり、リップル(揺らぎ)が無くいため、除去帯域近くの動きも滑らかです。※高次になれば傾斜は急になります
また、通過帯域は今から説明する他のフィルタと比較して線形な位相応答を示しており、周波数ごとの位相遅れの差もほとんど無いのが特徴です。

以下がサンプルコードとなり、フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力する関数です。
「ftype」数字によってLPF/HPF/BPFを選択できるようになっており、設定結果及び「fl」「fh」を任意の周波数/「order」次数を設定することでフィルタのパラメータを設定します。
関数「character_MA」がフィルタ特性の出力用関数
関数「character_Butter」が実際のフィルタ処理用関数になります。

dt = 0.005 #SamplingTime
fs = 1/dt  #周波数
ftype = 0 #lowpass=0, highpass=1, bandpass=2
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)#第四引数Trueならばanarog filter/ FalseならばDigital filter
        w, h = signal.freqz(b, a)#アナログフィルタはfreqs,デジタルフィルタはfreqz

        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)#第四引数Trueならばanarog filter/ FalseならばDigital filter
        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)  # MATLAB fvtool仕様
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)  # MATLAB fvtool仕様
plt.xlabel('Frequency [Hz]')
plt.ylabel('Group delay [samples]')
plt.grid()
plt.show()

実行結果は以下の通りです。
f:id:toku_dango:20210718133317p:plain
※サンプルコードはローパスフィルタ設定
グラフ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()

f:id:toku_dango:20210718133720p:plain
グラフ1つ目の緑プロットがButterworthフィルタ後の出力となります。※フィルタ処理前が赤プロット
グラフ2つ目がフィルタ後出力の周波数特性を示した時系列です。
グラフ3つ目の緑プロットが移動平均フィルタ後の周波数特性になります。
設計通り、20Hz付近から減衰するようなフィルタ特性を持っており、グラフ3つ目の緑プロットと赤プロットを比較して、設計通りのフィルタが作成できていることが分かります。

Chebyshevフィルタについて

通過帯域でのリップルが存在するが、Butterworseと比較して減衰頻度を大きく取れるフィルタ(第一種Chebyshevフィルタ)
減衰頻度をさらに大きく取れるように設計されたフィルタが第二種Chebyshevフィルタ

第一種Chebyshevフィルタ(フィルタ特性/フィルタ処理)

以下がサンプルコードとなり、フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力する関数です。
「ftype」数字によってLPF/HPF/BPFを選択できるようになっており、設定結果及び「fl」「fh」を任意の周波数/「order」次数を設定することでフィルタのパラメータを設定します。
第一種Chebyshevフィルタでは、「attenuation」でリップルについても調整可能です。

dt = 0.005 #SamplingTime
fs = 1/dt  #周波数
ftype = 0 #lowpass=0, highpass=1, bandpass=2
fl = None #カットオフ周波数(遮断周波数) 低周波数側
fh = 20 #カットオフ周波数(遮断周波数) 高周波数側
order = 4 #次数
attenuation = 5 #リップル(通過域最大損失量[dB])

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)#アナログフィルタはfreqs,デジタルフィルタはfreqz

        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)#第四引数Trueならばanarog filter/ FalseならばDigital filter

        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フィルタと比較して減衰が大きくなっているように見えます。
ただ、位相/群遅延プロットを見ると、周波数によって遅れが変わってきているので、使う際はこの辺を理解したうえで使う必要ありそうです。
f:id:toku_dango:20210718140925p:plain


フィルタ処理結果は以下の通りです。
f:id:toku_dango:20210718141244p:plain
ローパスフィルタとしてうまく機能はしてそうですが、リップルの影響か、sin波の強度が若干低下していますね。
この辺も理解したうえで使う必要ありそうです。

第二種Chebyshevフィルタ(フィルタ特性/フィルタ処理)

第一種Chebyshevフィルタよりも減衰頻度を大きく取れるように設計されたフィルタ
通過帯域にリップルがない代わりに、除去帯域で等リップル性を持つ。
以下がサンプルコードとなり、フィルタ評価に使用するサンプル信号とサンプル信号の周波数特性を出力する関数です。
「ftype」数字によってLPF/HPF/BPFを選択できるようになっており、設定結果及び「fl」「fh」を任意の周波数/「order」次数を設定することでフィルタのパラメータを設定します。
第二種Chebyshevフィルタでは、「attenuation」でリップルについても調整可能です。

dt = 0.005 #SamplingTime
fs = 1/dt  #周波数
ftype = 0 #lowpass=0, highpass=1, bandpass=2
fl = None #カットオフ周波数(遮断周波数) 低周波数側
fh = 20 #カットオフ周波数(遮断周波数) 高周波数側
order = 4 #次数
attenuation = 10 #minimum attenuation(阻止域最小減衰量[dB])

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)#アナログフィルタはfreqs,デジタルフィルタはfreqz

        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)

フィルタ特性は以下の通りです。
f:id:toku_dango:20210718144639p:plain

フィルタ処理結果は以下の通りです。
f:id:toku_dango:20210718142736p:plain


Besselフィルタについて

最後にBesselフィルタについて紹介します。
Butterworthと非常に似た特性を持つフィルタで、過渡応答性がより滑らかなフィルタ。
通過帯域では線形な位相特性を持つため、信号の遅れも周波数によらずある程度一定を保てることができます。
正直言うとButterworthとBesselフィルタとの使い分けはわかりませんでした。もう少し勉強します。

フィルタ特性/フィルタ処理
dt = 0.005 #SamplingTime
fs = 1/dt  #周波数
ftype = 0 #lowpass=0, highpass=1, bandpass=2
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)#アナログフィルタはfreqs,デジタルフィルタはfreqz

        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)#第四引数Trueならばanarog filter/ FalseならばDigital filter

        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)



フィルタ特性は以下の通りです。
f:id:toku_dango:20210718143350p:plain
確かにButterworthフィルタと比較して減衰域が滑らかになっているように見えますね。ただ、通過域と減衰域の境が滑らかになっており攻めたフィルタをしたときに境目付近の信号は減衰されてしまう可能性がありそうですね。


フィルタ処理結果は以下の通りです。
f:id:toku_dango:20210718143405p:plain

まとめ

今回はScipyを使ったディジタルフィルタの設計・評価で使えるコードを紹介しました。
正直位相特性も線形なButterworthが一番汎用性高いのではと思ってしまうのですが、使い分けはもう少し勉強します。