Scipyで多段のフィルタをかけ合わせる

背景

ScipyではMATLAB-likeにフィルタを設計するコマンド群が存在します。
計算結果は分母と分子の係数として出てくるのですが、これらを多段に重ねた結果を取得したい場合の対処法をメモします。

実装

numpyのpoly1dを使った多項式計算で簡単にフィルタの掛け合わせが実現できます。
公式の関数がぱっと見当たらなかったので作りましたが、無い訳は無い気がしているので見つけたら追記します。

なお、この方式だと極ゼロ相殺の判定ができないので自分で設計した覚えのないフィルタにこれを使用するのは推奨しません。

import scipy
import numpy as np

def series_filter(b_s,a_s):
    a = np.poly1d([1])
    b = np.poly1d([1])
    
    for ai, bi in zip(a_s,b_s):
        a *=  np.poly1d(ai)
        a *=  np.poly1d(bi)
    return b.coeffs, a.coeffs

使用例

例えばLowpassとノッチを組み合わせるなどすると下記のような離散フィルタを設計できます。

MATLABと同じfreqzというコマンドで上記の可視化ができるのでお試しあれ。

import matplotlib.pyplot as plt

b1,a1 = scipy.signal.butter(2,10, btype='low',analog=False,fs=1000)
b2,a2 = scipy.signal.butter(2,[20,25], btype='bandstop',analog=False,fs=1000)

b,a = series_filter([b1,b2],[a1,a2])

## 可視化
w, h = signal.freqz(b,a,fs=1000)
fig = plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.set_title('Digital filter frequency response')
ax1.plot(w, 20 * np.log10(abs(h)), 'b')
ax1.set_ylabel('Amplitude [dB]', color='b')
ax1.set_xlabel('Frequency [rad/sample]')
ax1.grid()
ax1.set_xscale('log')
ax1.set_xlim([0.1,100])
angles = np.unwrap(np.angle(h))/np.pi * 180.0
ax2=fig.add_subplot(2,1,2)
ax2.plot(w, angles, 'g')
ax2.set_ylabel('Angle (deg)', color='g')
ax2.grid()
ax2.axis('tight')
ax2.set_xscale('log')
ax2.set_xlim([0.1,100])
plt.show()