パラメトリックEQ

先日Juliaでフィルタの可視化をしたけれど、あの時点ではまだフィルタの設計はMatlabに頼っていた。今回はRobert Bristow-Johnson氏の名作をもとに、パラメトリックEQを作る。ただし、今日ここに紹介するのはピークフィルタのみ。

まず、先日のプログラムを関数にまとめておく。Juliaの流儀がいまいち分かっていないので、Python+SciPyだったりMatlabの雰囲気になっているかも。標本化周波数は44,100 Hzをデフォルトにしておいた。

function freqz(B, A, samprate=44100)
    # calculate frequency
    # (using logspace since we are going to use semilogx for plotting)
    fL = 10;
    fH = samprate/2;
    f = logspace(log10(fL), log10(fH), 1001);
    z = f / samprate * 2*pi;
    s = exp.(-1.0im * z);   # convert from digital frequency
    
    # calculate transfer function
    H = (B[1] + B[2].*s + B[3].*s.*s) ./ (A[1] + A[2].*s + A[3].*s.*s);
    
    subplot(2,1,1);
    semilogx(f, 20*log10.(abs.(H)));
    grid("on");
    ylabel("Gain (dB)");
    
    subplot(2,1,2);
    semilogx(f, angle.(H)*180/pi);
    grid("on");
    xlabel("Frequency (Hz)");
    ylabel("Phase (angle)");
end

Robert Bristow-Johnson氏の名作であるところの「Cookbook formulae for audio EQ biquad filter coefficients（の日本語訳）」から、Peaking EQの部分だけを関数にした。標本化周波数はここも44,100 Hz、Qは1/√2にしておいた。

function biquad_peak(fc, dBgain, Q=1/sqrt(2), samprate=44100)
    A = sqrt(10^(dBgain/20));
    w = 2 * pi * fc / samprate;
    s = sin(w);
    c = cos(w);
    Q = Q/A;   # adjustment suggested by RBJ (and matches with Zolzer's version) 2018-03-16
    a = s / (2 * Q);
    
    b0 =  1 + a*A;
    b1 = -2 * c;
    b2 =  1 - a*A;
    a0 =  1 + a/A;
    a1 = -2 * c;
    a2 =  1 - a/A;
    
    return([b0 b1 b2], [a0 a1 a2]);
end

で、Juliaから使うときには以下のようにする。

julia> using PyPlot
julia> B, A = biquad_peak(1000, 12, 2.0, 44100)
([1.14132 -1.97973 0.858678], [1.0355 -1.97973 0.964501])
julia> freqz(B, A);
julia> savefig("output.png");

すると以下のようなフィルタ形状がファイルで得られる。RBJ Cookbookをほぼそのまま変更なしで使えるのでありがたい。

f:id:amarui:20180318133846p:plain