DSP.jlに対応した低域フィルタ

引き続きJuliaネタです。前回は、Bristow-JohnsonのパラメトリックEQを自作のfreqz関数で可視化してみました。

marui.hatenablog.com

ただ、プログラミング言語を使う際には長いものには巻かれておいたほうが良いので、Juliaのエコシステムで言及されているJuliaDSPに対応したものに更新したいと思います。フィルタ設計などで特に使うのはDSP.jlです。

まず必要なパッケージを読み込みます。音の読み書きにLibSndFileと、可視化のためにPlotsも使います。

using DSP, LibSndFile, SampledSignals, Plots
pyplot();

今回はBristow-Johnsonのレシピから低域フィルタを拝借します。これまではフィルタ係数を配列に入れて戻していましたが、それをBiquad型にするだけでDSP.jlと相性が良くなります。Biquad型はa0=1.0を想定しているので、各係数をa0で割っています。

function biquad_lpf(fc, Q=1/sqrt(2), samprate=44100)
    w = 2 * pi * fc / samprate;
    s = sin(w);
    c = cos(w);
    a = s / (2 * Q);

    b0 = (1 - c) / 2;
    b1 =  1 - c;
    b2 = (1 - c) / 2;
    a0 =  1 + a;
    a1 = -2 * c;
    a2 =  1 - a;

    return(Biquad(b0/a0, b1/a0, b2/a0, a1/a0, a2/a0));
end

DSP.jlのfreqzを使ってフィルタの応答を計算してグラフ描画します。

fs = 48000;
fltr = biquad_lpf(1000, 2, fs);
f = 20:20000;
A = freqz(fltr, f, fs);
res_power = plot(f, 20*log10.(abs.(A)),
    xscale=:log10, ylim=(-45, 15), legend=false,
    xticks=([31.5, 63, 125, 250, 500, 1000, 2000, 4000, 8000, 16000],
        ["31.5", "63", "125", "250", "500", "1k", "2k", "4k", "8k", "16k"]),
    xlabel="Frequency (Hz)", ylabel="Gain (dB)");
res_phase = plot(f, angle.(A),
    xscale=:log10, ylim=(-pi, +pi), legend=false,
    xticks=([31.5, 63, 125, 250, 500, 1000, 2000, 4000, 8000, 16000],
        ["31.5", "63", "125", "250", "500", "1k", "2k", "4k", "8k", "16k"]),
    xlabel="Frequency (Hz)", ylabel="Phase (rad)");
plot(res_power, res_phase, layout=grid(2,1))

音ファイルへのフィルタの適用は以下のようにやります。音ファイルを読み込み、filt関数でフィルタをかけてから保存しています。

snd = load("input.wav");
y = filt(fltr, snd.data);
buf = SampleBuf(y, snd.samplerate);
save("output.wav", buf);

クラシックギター音に適用してみました。

【適用前】

【適用後】

ちゃんとフィルタがかかっているように見えます。