【scipy】でSOS形式のIIRフィルターの周波数応答を見たい

IIRフィルターについてはscipyのドキュメントが詳しいのでそちらを一読がおすすめ

■ 参考 : Signal Processing with SciPy: Linear Filters(Warren Weckesser)

IIRフィルターには3つの表し方の形式がある
・伝達関数
・zero, pole, gain(ZPK方式)
・second order sections (SOS形式)

最近までSOS形式のことをSOSフィルターだと思ってたけど、これはフィルターではないっぽい

使うのはscipyのsosfreqs
これで、SOS形式の周波数応答が書ける

from scipy.signal import butter, sosfreqz
import numpy as np
import matplotlib.pyplot as plt

fs=16384
fs_nyq = fs*0.5
low = 100 / fs_nyq
high = 500 /fs_nyq

sos = butter(20, [low, high], btype="bandpass", output='sos')
w, h = sosfreqz(sos, worN=8000)
plt.plot((fs_nyq/np.pi)*w, np.abs(h))
plt.plot((fs_nyq/np.pi)*w, np.angle(h))
ax = plt.gca()
ax.set_xlabel("frequency [Hz]")
ax.set_xlim([0, 1000])
plt.show()

fs=16384
fs_nyq = fs*0.5
low = 100 / fs_nyq
high = 500 /fs_nyq

sos = butter(5, [low, high], btype="bandpass", output='sos')
w, h = sosfreqz(sos, worN=8000)
plt.plot((fs_nyq/np.pi)*w, np.abs(h))
plt.plot((fs_nyq/np.pi)*w, np.angle(h))
ax = plt.gca()
ax.set_xlabel("frequency [Hz]")
ax.set_xlim([0, 1000])
plt.show()

この例は100Hz ~ 500Hzのバンドパスフィルター
次数は20次と5次

■ 参考 : scipy.signal.butter

■ 参考 : scipy.signal.ellip

SOSフィルターの係数がわかっている場合

from scipy.signal import butter, sosfreqz, ellip
import numpy as np
import matplotlib.pyplot as plt
import math

import matplotlib.ticker as ptick
from matplotlib.ticker import ScalarFormatter, NullFormatter


import matplotlib as mpl

mpl.rcParams['xtick.labelsize'] = 20
mpl.rcParams['ytick.labelsize'] = 20
mpl.rcParams["axes.labelsize"] = 20
mpl.rcParams["axes.titlesize"] = 20
mpl.rcParams["legend.fontsize"] = 20

# filter coefficient
# order : [b0, b1, b2, a0, a1, a2]
#
label=band pass 0.03-0.1 Hz"
freq_range = [0.02, 0.2]
coeff1 = [1, -1.999998669372317, +0.999999999999999, 1, -1.999946242188092, +0.999946281201766]
coeff2 = [1, -1.999999767511141, +1.000000000000011, 1, -1.999955898604725, +0.999955963180191]
coeff3 = [1, -1.999999999400776, +0.999999999999995, 1, -1.999961099937338, +0.999961120373777]
coeff4 = [1, -1.999999853828710, +0.999999999999960, 1, -1.999976229681027, +0.999976314187963]
coeff5 = [1, -1.999999996570461, +1.000000000000034, 1, -1.999980731688244, +0.999980744035228]
coeff6 = [1, -1.999999994545118, +0.999999999999935, 1, -1.999992076212929, +0.999992085647952]
coeff7 = [1, -1.999999873305222, +1.000000000000034, 1, -1.999992919697037, +0.999993013663749]
coeff8 = [1, -1.999999993706660, +1.000000000000031, 1, -1.999997892111085, +0.999997900596364]

coeff = [coeff1, coeff2, coeff3, coeff4,
coeff5, coeff6, coeff7, coeff8]
sos = np.array(coeff, dtype="float64")
print(coeff)

#
# filter response
#
w, h2 = sosfreqz(sos, worN=1000000) # worN decided the frequency resolution

#
# plot
#

plt.figure(figsize=[12, 6])
plt.plot((fs_nyq/np.pi)*w / factor, np.abs(h), label="hoge", color="magenta")
ax = plt.gca()
ax.set_xlabel("frequency [Hz]")
ax.set_ylabel("frequency response")
ax.set_xlim([1e-4, fs_nyq])
ax.set_title("with correction of down-sampling")
plt.xscale("log")
plt.yscale("log")
plt.legend(loc="lower right")
plt.gca().xaxis.set_major_formatter(plt.FormatStrFormatter('%.2f'))#y軸小数点以下3桁表示
plt.show()