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を使って高速に計算するのかの違いである。