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);
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
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
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
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
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
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
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
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
この項はUdo Z lzer “Digital Audio Signal Processing” (1997)の第5章から取ってきたフィルタ係数。Robert Bristow-Johnsonの係数と同じ方法で使用できる。
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
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
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
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
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
入力音を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