このエントリはJulia Advent Calendar 2020の25日目のものです。Juliaを音に関する研究によく使っているので、その中から一つ紹介したいと思います。
音の加工をするときにフィルタというものを使うことがあります。フィルタは、特定の周波数の音を除去したり強調したりするエフェクタです。フィルタを使って音を周波数帯域ごとに分割する、フィルターバンクというものもあります。周波数帯域ごとにどんな音が含まれているのかを調べるときに使えますし、一部の帯域だけを強調したりするグラフィック・イコライザとして使用することもできます。
グライコとして使うのであれば、フィルタバンクで分割した音をすべて再合成すると元の音に戻ることが必須になります。元に戻るということを必須条件にしているフィルタバンクには直交ミラーフィルタ（QMF）というものもありますが、今回はもっと簡単な原理のものを紹介します。
Raised-cosine
Raised-cosineという、余弦波（コサイン波）の一部を使うものがあります。下の図には灰色の線で余弦波を描いています。raised-cosineはその余弦波のうち振幅が最も低くなったところからの1周期（赤線部分）を抜き出して、振幅を[0, 1]の範囲にスケール（raise）したものです。
これは音のフェードインやフェードアウトのときにもスムーズなカーブとして使うことができます。とくにクロスフェード（音Aをフェードアウトさせながら音Bをフェードインさせるなど）のときには、次に説明する性質が便利に利用できます。
下の図のように、フェードアウトに入っていくタイミングとフェードインのタイミングが合っているときには、位相がπ/2ずれます。すると、クロスフェード部分はcosθとcos(θ - π/2)（つまりsinθ）の波となっていることが分かります。この図の縦軸は振幅ですが、エネルギ（振幅の2乗に比例）で考えたほうが、聴感では自然に感じられます。エネルギは振幅の2乗に比例するので、cosとsinをそれぞれ2乗して足し合わせるとエネルギの計算になります。ここで cos²θ + sin²θ = 1 という性質から、クロスフェードのどのタイミングでもエネルギは一定のままになることが分かります。つまり、音の大きさがあるていど等しく聞こえることが期待できます。*1
Raised-cosineを使ったフィルタバンク
このRaised-cosineをフィルタバンクに応用します。下のグラフのように、各バンドのフィルタ形状がraised-cosineのかたちになっているものだと、前述のように隣り合うバンドを加算したときに中間周波数部分のエネルギが一定になります。つまり帯域ごとに分割した音を加算すると、また元の信号に戻せるようなフィルターバンクというわけです。このフィルタ形状はなんらかの規格にそったものではないので、他の方法と比べるときには不都合があるかもしれませんが、グラフィック・イコライザなど様々な用途に使えそうです。
フィルタバンクを実装するときには、(1)入力音をフーリエ変換し、(2)raised-cosine窓をかけて、(3)個別のバンドを逆フーリエ変換で音に戻す、というのが全体の流れです。
高速フーリエ変換した信号は線形周波数で振幅・位相が並びます（このあたりは以前のエントリを参照して下さい）。ここで面倒なのが、raised-cosineフィルタバンクをグラフィック・イコライザのように使用したいとすると、周波数を対数周波数にしないといけないということです。たとえば、1 kHzを中心周波数にしたバンドパスフィルタを見ると、低い周波数側は500〜1000 Hzの500 Hz幅なのに対して、高い周波数側は1000〜2000 Hzと1000 Hz幅になります。対数周波数上ではフィルタ形状が左右対称になるのですが、線形周波数で考えると高域側へ尻尾が伸びたかたちになるのです。
1000 Hzを中心周波数とした1オクターブ幅のバンドパスフィルタは、500〜2000 Hzをフィルタの両端にします。その計算は以下のようにします。中心周波数をfc、1/Nオクターブ幅のフィルタの両端の周波数（fLとfH）は次の式で計算できます。*2
fL = fc * 2^(-1 / octave); fH = fc * 2^(+1 / octave);
次に、窓Wを準備します。通過帯域にはraised-cosineを入れ、それ以外をゼロにしたいです。そこで、まず窓全体をゼロで埋めておきます。
W = zeros(round(Int, length(xx) / 2));
通過帯域はfL〜fHなので、その範囲を調べます。以下の式ではfにDC〜ナイキスト周波数を入れ、そのうち通過帯域に入っている部分をtrue、それ以外をfalseにしたffを作ります。
f = range(0, stop = samprate / 2, length = length(W)); ff = (fl .<= f) .& (f .<= fh);
ccにraised-cosineを作ります。ここでlog2を使っているのは、前述の対数周波数軸で左右対称にするためです。Wの通過帯域部分に入れて窓が完成です。
cc = (cos.(log2.(f[ff] / fc) * pi * octave) .+ 1) / 2; W[ff] = cc;
さて、最も低い周波数帯域と最も高い周波数帯域（上記フィルタバンク図で31.5 Hz帯と16 kHz帯）をバンドパスフィルタにすると、それより低い・高いところにある両端の周波数帯域の信号が消えてしまいます。そこで、最低・最高周波数帯域は以下のように別の対応をします。（ここでnはバンド番号、各中心周波数以下・以上を1にしています）
if n == 1 W[f .< fc] .= 1; elseif n == length(cfreqs) W[fc .< f] .= 1; end
逆フーリエ変換をするときには、正の周波数と負の周波数の両側を考慮する必要があります。念のため、DC成分とナイキスト周波数成分はゼロにします。負の周波数側は複素共役になるはずですが、今回の窓の要素はすべて 0+0im か 1+0im なのでconjしなくてOKです。
W = [W; reverse(W)]; W[1] = 0; W[round(Int, length(W) / 2) + 1] = 0;
まとめたもの
長々と説明しましたが、以下が関数全体です。
""" 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 = [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/N octave 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
以下のようにすると、hoge.wavに入っている音データを1/3オクターブごとに分割して、hoge2.wavに保存できます。hoge2.wavには30チャンネルの音が入っていますが、各チャンネルはcfreqsの中心周波数で1/3オクターブ分割された音です。filterbank(x, fs, 12)とすれば1/12オクターブ（＝半音）単位での分析もできます。ただ、全体の基準にしている周波数は1000 Hzなので、音楽分析の用途で使うには440 Hzを基準にするなどのプログラム修正が必要になるかと思います。
using FileIO: load, save import LibSndFile using SampledSignals snd = load("hoge.wav"); x = snd.data; fs = snd.samplerate; out, cfreqs = filterbank(x, fs, 3); snd2 = SampleBuf(out, fs); save("hoge2.wav", snd2);
