MatlabのfreqzもどきをJuliaで

MatlabやOctaveにはfreqzという関数があり、FIR/IIRフィルタの周波数特性・位相特性を計算しグラフ表示してくれます。PythonでもSciPy.signalに同名の関数があります。似たことをJuliaでやってみました。

今回は双二次フィルタを見てみましょう。標本化周波数44,100 Hz、中心周波数1,000 Hz、ゲイン12 dB、Q=2.0のピークフィルタを準備しました。

samprate = 44100;
B = [1.1417, -1.9797, 0.8583];
A = [1.0356, -1.9797, 0.9644];

特性を見たい周波数範囲を設定し、それをzさらにsへと変換して使います。ここでは16 Hz〜ナイキスト周波数を256段階の対数分割してみました。

# calculate frequency
# (using logspace since we are going to use semilogx for plotting)
fL = 16;
fH = samprate/2;
f = logspace(log10(fL), log10(fH), 256);
z = f / sampfreq * 2*pi;
s = exp.(-1.0im * z);   # convert from digital frequency

次に、各周波数での伝達関数の値を計算します。下の式に代入して計算するだけ。

[tex: H=\frac{b_0 + b_1 s + b_2 s²}{a_0 + a_1 s + a_2 s²}]

# calculate transfer function
H = (B[1] + B[2].*s + B[3].*s.*s) ./ (A[1] + A[2].*s + A[3].*s.*s);

Polynomialパッケージを使って汎用性を高めても良いのですが、僕の考えている使用用途は双二次フィルタだけなので、このままで。

さて、振幅特性をグラフにしてみましょう。abs.()で絶対値を計算してからデシベルにします。

using PyPlot
semilogx(f, 20*log10.(abs.(H)));
grid();
xlabel("Frequency (Hz)");
ylabel("Gain (dB)");
savefig("freqz_abs.png")

次に位相特性もグラフにしてみましょう。今度はangle.()で位相を抜き出し、ラジアンを角度に変換しています。

semilogx(f, angle.(H)*180/pi);
grid();
xlabel("Frequency (Hz)");
ylabel("Phase (angle)");
savefig("freqz_phase.png")

Matlabのfreqzのグラフと見比べてみましたが、だいたい同じ結果になっているようでした。Matlab版は横軸が線形になっているので見た目が違いますが。