scipy.signal.freqz相当の処理を自分で書く

これまで仕事でデジタルフィルタの周波数特性を見るためにscipy.signal.freqzを使っていたのだが、連続して実行するとどうにも遅い。 ipythonでは%runで外部スクリプトを実行する際に-pを付けるとプロファイル結果を表示することが可能であるが、freqzを使っているスクリプトで簡単に測ってみたところ、どうやらfreqzの中で使われているnumpy.polyvalが重いことが分かった。そこで、デジタルフィルタの周波数特性の計算をfftで行うことを検討した。

結論から書くと、係数ベクトルが1次元の時は以下のような関数で代替可能である。正規化周波数の終端はscipy.signal.freqzと同じくπまでとしているため、wを周波数に直したい場合は(FS / (2.0 * np.pi))を掛ける必要がある。

import numpy as np
import scipy as sp
import scipy.fftpack as fft

def my_freqz(b, a=[1], worN=None):
    lastpoint = np.pi
    N = 512 if worN is None else worN
    w = np.linspace(0.0, lastpoint, N, endpoint=False)
    h = fft.fft(b, 2 * N)[:N] / fft.fft(a, 2 * N)[:N]
    return w, h

b = np.random.randn(1024)
a = np.random.randn(1024)

w, h = sg.freqz(b, a, 512)
my_w, my_h = my_freqz(b, a, 512)

np.allclose(h, my_h)  # -> True
np.sum(h - my_h) # -> (3.6308811859897538e-12-1.5773658421092129e-11j)

計算方法が違うため、誤差はあるものの非常に小さいため十分代替可能である。 ポイントはfftを2N点で計算して前半N点だけを使うところか。なおwの計算方法はscipy.signal.freqzと同じである。

折角なのでipython上で%timeitを使って時間を計測してみると、

%timeit -n1000 sg.freqz(b, a, 512)
1000 loops, best of 3: 18.7 ms per loop

%timeit -n1000 my_freqz(b, a, 512)
1000 loops, best of 3: 72.6 µs per loop

となり、速度的には上記実装の方が圧倒的に速い。ただし、scypy.signal.freqzは係数に対してatleast_1dを実行したり引数のチェックなども入っているため、実際には完全に同じ条件下での比較ではない。

なお、numpy.polyvalを使った計算とfftを使った計算結果が同じになる理由については、デジタルフィルタの周波数特性の定義式を素直にDFTとして計算するのか、FFTを使って高速に計算するのかの違いである。