Subsections

2 エフェクタ

2.1 イコライザ(Robert Bristow-Johnson版)

Robert Bristow-Johnsonの“Audio EQ Cookbook”の双二次フィルタをそのままコードにした。次のようにして使うことができる。(入力ファイルが1chと仮定している)

using DSP
using SampledSignals
using FileIO: load, save
import LibSndFile
snd = load("test.wav");
fltr = biqual_lpf_rbj(1000, 2, snd.samplerate);
y = filt(fltr, snd.data);
snd2 = SampleBuf(y, snd.samplerate);
save("test2.wav", snd2);

2.1.1 低域フィルタ(LPF)

function biquad_lpf_rbj(fc, Q = 1 / sqrt(2), samprate = 48000)
  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

2.1.2 高域フィルタ(HPF)

function biquad_hpf_rbj(fc, Q = 1 / sqrt(2), samprate = 48000)
  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

2.1.3 帯域フィルタ(BPF)

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

  b0 =  a;
  b1 =  0;
  b2 = -a;
  a0 =  1 + a;
  a1 = -2 * c;
  a2 =  1 - a;

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

2.1.4 ノッチフィルタ

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

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

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

2.1.5 オールパスフィルタ

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

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

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

2.1.6 ピークフィルタ

function biquad_peak_rbj(fc, dBgain, Q = 1 / sqrt(2), samprate = 48000)
  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)
  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 Biquad(b0/a0, b1/a0, b2/a0, a1/a0, a2/a0);
end

2.1.7 ローシェルフフィルタ

function biquad_lowshelf_rbj(fc, gain, S = 1, samprate = 48000)
  A = sqrt(10^(gain / 20));
  w = 2 * pi * fc / samprate;
  s = sin(w);
  c = cos(w);
  a = s / 2 * sqrt((A + 1 / A) * (1 / S - 1) + 2);
  b = 2 * sqrt(A) * a;

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

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

2.1.8 ハイシェルフフィルタ

function biquad_highshelf_rbj(fc, gain, S = 1, samprate = 48000)
  A = sqrt(10^(gain / 20));
  w = 2 * pi * fc / samprate;
  s = sin(w);
  c = cos(w);
  a = s / 2 * sqrt((A + 1 / A) * (1 / S - 1) + 2);
  b = 2 * sqrt(A) * a;

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

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

2.2 イコライザ(Udo Zlzer版)

この項はUdo Zlzer “Digital Audio Signal Processing” (1997)の第5章から取ってきたフィルタ係数。Robert Bristow-Johnsonの係数と同じ方法で使用できる。

2.2.1 低域フィルタ

function biquad_lpf_zolzer(fc, samprate = 48000)
  K = tan(pi * fc / samprate);
  b0 =   K^2 / (1 + sqrt(2) * K + K^2);
  b1 = 2 * K^2 / (1 + sqrt(2) * K + K^2);
  b2 =   K^2 / (1 + sqrt(2) * K + K^2);
  a0 = 1;
  a1 = 2 * (K^2 - 1) / (1 + sqrt(2) * K + K^2);
  a2 = (1 - sqrt(2) * K + K^2) / (1 + sqrt(2) * K + K^2);

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

2.2.2 高域フィルタ

function biquad_hpf_zolzer(fc, samprate = 48000)
  K = tan(pi * fc / samprate);
  b0 =  1 / (1 + sqrt(2) * K + K^2);
  b1 = -2 / (1 + sqrt(2) * K + K^2);
  b2 =  1 / (1 + sqrt(2) * K + K^2);
  a0 = 1;
  a1 = 2 * (K^2 - 1) / (1 + sqrt(2) * K + K^2);
  a2 = (1 - sqrt(2) * K + K^2) / (1 + sqrt(2) * K + K^2);

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

2.2.3 ローシェルフフィルタ

function biquad_lowshelf_zolzer(fc, gain, samprate = 48000)
  K = tan(pi * fc / samprate);
  if gain >= 0.0
    V0 = 10^(gain / 20);
    b0 = (1 + sqrt(2 * V0) * K + V0 * K^2) / (1 + sqrt(2) * K + K^2);
    b1 =                2 * (V0 * K^2 - 1) / (1 + sqrt(2) * K + K^2);
    b2 = (1 - sqrt(2 * V0) * K + V0 * K^2) / (1 + sqrt(2) * K + K^2);
    a0 = 1;
    a1 =                     2 * (K^2 - 1) / (1 + sqrt(2) * K + K^2);
    a2 =           (1 - sqrt(2) * K + K^2) / (1 + sqrt(2) * K + K^2);
  else
    V0 = 10^(-gain / 20);
    b0 =           (1 + sqrt(2) * K + K^2) / (1 + sqrt(2 * V0) * K + V0 * K^2);
    b1 =                     2 * (K^2 - 1) / (1 + sqrt(2 * V0) * K + V0 * K^2);
    b2 =           (1 - sqrt(2) * K + K^2) / (1 + sqrt(2 * V0) * K + V0 * K^2);
    a0 = 1;
    a1 =                2 * (V0 * K^2 - 1) / (1 + sqrt(2 * V0) * K + V0 * K^2);
    a2 = (1 - sqrt(2 * V0) * K + V0 * K^2) / (1 + sqrt(2 * V0) * K + V0 * K^2);
  end

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

2.2.4 ハイシェルフフィルタ

function biquad_highshelf_zolzer(fc, gain, samprate = 48000)
  K = tan(pi * fc / samprate);
  if gain >= 0.0
    V0 = 10^(gain / 20);
    b0 =     (V0 + sqrt(2 * V0) * K + K^2) / (1 + sqrt(2) * K + K^2);
    b1 =                    2 * (K^2 - V0) / (1 + sqrt(2) * K + K^2);
    b2 =     (V0 - sqrt(2 * V0) * K + K^2) / (1 + sqrt(2) * K + K^2);
    a0 = 1;
    a1 =                     2 * (K^2 - 1) / (1 + sqrt(2) * K + K^2);
    a2 =           (1 - sqrt(2) * K + K^2) / (1 + sqrt(2) * K + K^2);
  else
    V0 = 10^(-gain / 20);
    b0 =           (1 + sqrt(2) * K + K^2) / (V0 + sqrt(2 * V0) * K + K^2);
    b1 =                     2 * (K^2 - 1) / (V0 + sqrt(2 * V0) * K + K^2);
    b2 =           (1 - sqrt(2) * K + K^2) / (V0 + sqrt(2 * V0) * K + K^2);
    a0 = 1;
    a1 =                2 * (K^2 / V0 - 1) / (1 + sqrt(2 / V0) * K + K^2 / V0);
    a2 = (1 - sqrt(2 / V0) * K + K^2 / V0) / (1 + sqrt(2 / V0) * K + K^2 / V0);
  end

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

2.2.5 ピークフィルタ

function biquad_peak_zolzer(fc, gain, Q = 1 / sqrt(2), samprate = 48000)
  K = tan(pi * fc / samprate);
  if gain >= 0.0
    V0 = 10^(gain / 20);
    b0 = (1 + V0 / Q * K + K^2) / (1 + 1 / Q * K + K^2);
    b1 =          2 * (K^2 - 1) / (1 + 1 / Q * K + K^2);
    b2 = (1 - V0 / Q * K + K^2) / (1 + 1 / Q * K + K^2);
    a0 = 1;
    a1 =          2 * (K^2 - 1) / (1 + 1 / Q * K + K^2);
    a2 =  (1 - 1 / Q * K + K^2) / (1 + 1 / Q * K + K^2);
  else
    V0 = 10^(-gain / 20);
    b0 =  (1 + 1 / Q * K + K^2) / (1 + V0 / Q * K + K^2);
    b1 =          2 * (K^2 - 1) / (1 + V0 / Q * K + K^2);
    b2 =  (1 - 1 / Q * K + K^2) / (1 + V0 / Q * K + K^2);
    a0 = 1;
    a1 =          2 * (K^2 - 1) / (1 + V0 / Q * K + K^2);
    a2 = (1 - V0 / Q * K + K^2) / (1 + V0 / Q * K + K^2);
  end

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

2.3 フィルタバンク

入力音をFFTして周波数領域にもっていき、raised-cosine窓を使って帯域分割するものです。raised-cosineが相補的なので、帯域分割したデータを加算すると元データに戻ります。元データに戻すために、最低帯域と最高帯域についてはカットオフ周波数以下・以上のすべてを入れてあります。たとえば、中心周波数31.5 Hzが最低の1オクターブ帯域であれば、上側カットオフ周波数はだいたい44.5 Hzになるので、44.5 Hz以下の音は全てその帯域に含むことにしています。

filterbank(【音データのベクトル】, 【サンプリング周波数】, 【1/オクターブ幅】)として使用します。【1/オクターブ幅】なので、1を指定したら1-oct、3を指定したら1/3-octになります。

using FFTW
using SampledSignals
using FileIO: load, save
import LibSndFile
snd = load("test.wav");
(y, freqs) = filterbank(snd.data, snd.samplerate, 1);
snd2 = SampleBuf(y, snd.samplerate);
save("test2.wav", snd2);   # 分割した帯域ぶんのチャンネル数になります

コードは以下です。

"""
    filterbank(x, samprate[, N=1])

Filter input signal `x` using 1/N-octave filterbank

`out, cfreqs = filterbank(x, samprate)` create filtered signals `out`
from `x`. Columns of `out` is a filtered signal with corresponding
center frequency in `cfreqs`. By default, 1/1-octave band filters are
used. Filtered signals at the lowest frequency is processed with
low-pass filter, the highest frequency is processed with high-pass
filter, and all frequencies in between are processed with band-pass
filters of raised cosine frequency window. This yields to (almost)
perfect reconstruction of the original signal, such that 
“`
out, cfreqs = filterbank(x, samprate)
z = sum(out, 1)
“`
creates `z` almost the same as `x`.

`out, cfreqs = filterbank(x, samprate, N)` creates and filters using
1/N-octave band filters.
"""
function filterbank(x, samprate = 48000, octave=1)
  # prepare center frequencies
  if octave == 1
    cfreqs = [31.5 63 125 250 500 1000 2000 4000 8000 16000 31500 \
              63000 125000 250000 500000 1000000];
  elseif octave == 3
    cfreqs = [16 20 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 \
              800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 \
              10000 12500 16000 20000 25000 31500 40000 50000 63000 \
              80000 100000 125000 160000 200000 250000 315000 400000 \
              500000 630000 800000 1000000];
  else
    cfreqs = 1000 * 2 .^ ((-10 * octave:10 * octave) / octave);
  end
  cfreqs = cfreqs[(20 .<= cfreqs) .& (cfreqs .<= samprate / 2)];

  # create 1/oct-band bandpass filter
  if length(x) % 2 == 1
    x = [x ; 0.0];
  end
  out = zeros((length(x), length(cfreqs)));
  xx = fft(x);
  for n in 1:length(cfreqs)
    fc = cfreqs[n];
    fl = fc * 2^(-1 / octave);
    fh = fc * 2^(+1 / octave);
    W = zeros(round(Int, length(xx) / 2));
    f = range(0, stop = samprate / 2, length = length(W));
    ff = (fl .<= f) .& (f .<= fh);
    cc = (cos.(log2.(f[ff] / fc) * pi * octave) .+ 1) / 2;
    W[ff] = cc;
    if n == 1
      W[f .< fc] .= 1;
    elseif n == length(cfreqs)
      W[fc .< f] .= 1;
    end
    W = [W; reverse(W)];
    W[1] = 0;
    W[round(Int, length(W) / 2) + 1] = 0;

    # apply filter
    y = real(ifft(xx .* W));
    out[:, n] = y;
  end
  return(out, cfreqs);
end



MARUI Atsushi
2022-04-20