04.フィルタ

4.1 フィルタの料理法

音を作る上でとっても重要なフィルタだが、デジタルフィルタという奴はあまりに自由度があってどう設計するべきか悩んでしまうのだ。が、デジタルフィルタの中でBiQuad型という奴は割合特性をコントロールしやすく、割合自由が効くという丁度いい感じなので色々と固まったノウハウがあるのだ。このBiQuadフィルタの係数の設計について Robert Bristow-Johnson という人が書いた有名な”Cookbook formulae for audio EQ biquad filter coefficients,”という文書があって、RBJ cookbookとか呼ばれている。

とっても役に立つので意訳気味に翻訳した。
原典はhttp://www.musicdsp.org/files/Audio-EQ-Cookbook.txtを参照してくれ。ちなみにデジタルフィルタの基本的な部分についてはhttp://www.digitalfilter.com/jpindex.htmlこのあたりが参考になるぞ。

Biquadフィルタの料理法

by Robert Bristow-Johnson

色んなフィルタのタイプについて解説するけど、どのフィルタの伝達関数もアナログな原型をバイリニア変換(BLT)してデジタル化したもんだぜ！このバイリニアの周波数変換は周波数の補正（これはBLTを使う時に必要ないわゆる”prewarp”という奴）とバンド幅の再調整（バンド幅はBLTでアナログからデジタルに変換した時に圧縮されるからな）を考慮してるぜ！

まず最初に、Bi-Quad伝達関数の基本の定義だ！

         b0 + b1*z^-1 + b2*z^-2
H(z) = ------------------------    (Eq 1)
        a0 + a1*z^-1 + a2*z^-2

この式では５個ではなく６個の係数を使うぜ！　アーキテクチャによっては、a0が１になるようにノーマライズしてb0も（多分）1になるぞ。（全体のゲインによるがな！）。すると伝達関数は次のような感じだ！

        (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2
H(z) = ---------------------------------------       (Eq 2)
           1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

または

                  1 + (b1/b0)*z^-1 + (b2/b0)*z^-2
H(z) = (b0/a0) * ---------------------------------      (Eq 3)
                  1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

一番直接的に実装するなら、Eq2の式がいいぜ！
つまりこういう事だ！！ y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2] (Eq 4) （注：　ここでx[n]は入力で、x[n-1],x[n-2]が入力をディレイさせたもの、y[n]が出力で、y[n-1],y[n-2]は出力をディレイさせたもの、という事だな。　この式をそのままプログラムにすればいいんだぜ！！！）

ＤＳＰの56Kという奴とか、似たような奴とか、浮動小数点を使う時には、これがベストで簡単だぜ！
じゃあ次に、ユーザー（あんた）が決めなくちゃいけないパラメータについてだ！

Fs （サンプリング周波数）
f0 （中心周波数またはカットオフ周波数またはシェルフの中間周波数。フィルタータイプによる。フィルターが”注目”する周波数）
dBGain （ピーキング、シェルビングの時にだけ使用）
Q （電気工学でいうフィルターのQだ（ただし、ピーキングEQの場合はA*Qで電気工学的なQと同じだ！これは、NデシベルカットしてからNデシベルブーストしたらフラットに戻るようにしてるためだよ））
Qの代わりにBWでも可。バンド幅（単位はオクターブ）だ（BPF/ノッチなら-3dB落ちるところ。ピーキングEQならdBGainが半分になる中間点の所の幅だ）
更にQの代わりのS。シェルビングEQの場合に使う”スロープ”のパラメータだ！　もしS=1ならスロープは可能な限り急峻な単調増加／減少になるぞ。スロープはdB/Octで表すと、f0/FsとdBgainが同じならSに比例するぞ。

よし、じゃあ、途中で使う変数を幾つか計算をするぜ！

A = sqrt( 10^(dBgain/20) )=10^(dBgain/40) (ピーキングとシェルビングEQの時だけな)
w0 = 2*pi*f0/Fs
cos(w0)
sin(w0)
alpha = sin(w0)/(2*Q) (Q を使う時) = sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) ) (BW を使う時) = sin(w0)/2 * sqrt( (A + 1/A)*(1/S – 1) + 2 ) (S を使う時)

　　　　バンド幅とＱの関係は、　 1/Q = 2*sinh(ln(2)/2*BW*w0/sin(w0)) (BLTのデジタルフィルタ) または 1/Q = 2*sinh(ln(2)/2*BW) (原型のアナログフィルタ) 　　　　シェルフのスロープとＱの関係は 1/Q = sqrt((A + 1/A)*(1/S – 1) + 2) 2*sqrt(A)*alpha = sin(w0) * sqrt( (A^2 + 1)*(1/S – 1) + 2*A ) 　　　　この式はシェルビングＥＱの時に便利な中間変数だ！で、最後に必要なフィルタータイプにあわせた係数の計算だ！（正規化された周波数に対するアナログの原型、H(s)、と対応する係数だ！）

LPF:  H(s) = 1 / (s^2 + s/Q + 1)
            b0 =  (1 - cos(w0))/2
            b1 =   1 - cos(w0)
            b2 =  (1 - cos(w0))/2
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

HPF:  H(s) = s^2 / (s^2 + s/Q + 1)
            b0 =  (1 + cos(w0))/2
            b1 = -(1 + cos(w0))
            b2 =  (1 + cos(w0))/2
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

BPF: H(s) = s / (s^2 + s/Q + 1)
         (constant skirt gain, peak gain = Q)

            b0 =   sin(w0)/2  =   Q*alpha
            b1 =   0
            b2 =  -sin(w0)/2  =  -Q*alpha
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

BPF:  H(s) = (s/Q) / (s^2 + s/Q + 1) 
     (constant 0 dB peak gain)

            b0 =   alpha
            b1 =   0
            b2 =  -alpha
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

notch:  H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
            b0 =   1
            b1 =  -2*cos(w0)
            b2 =   1
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

APF:  H(s) = (s^2 - s/Q + 1) / (s^2 + s/Q + 1)

            b0 =   1 - alpha
            b1 =  -2*cos(w0)
            b2 =   1 + alpha
            a0 =   1 + alpha
            a1 =  -2*cos(w0)
            a2 =   1 - alpha

peakingEQ:  H(s) = (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1)

            b0 =   1 + alpha*A
            b1 =  -2*cos(w0)
            b2 =   1 - alpha*A
            a0 =   1 + alpha/A
            a1 =  -2*cos(w0)
            a2 =   1 - alpha/A

lowShelf:
  H(s) = A * (s^2 + (sqrt(A)/Q)*s + A)/(A*s^2 + (sqrt(A)/Q)*s + 1)

            b0 =    A*( (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha )
            b1 =  2*A*( (A-1) - (A+1)*cos(w0)                   )
            b2 =    A*( (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha )
            a0 =        (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha
            a1 =   -2*( (A-1) + (A+1)*cos(w0)                   )
            a2 =        (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha

highShelf:
  H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A)

            b0 =    A*( (A+1) + (A-1)*cos(w0) + 2*sqrt(A)*alpha )
            b1 = -2*A*( (A-1) + (A+1)*cos(w0)                   )
            b2 =    A*( (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha )
            a0 =        (A+1) - (A-1)*cos(w0) + 2*sqrt(A)*alpha
            a1 =    2*( (A-1) - (A+1)*cos(w0)                   )
            a2 =        (A+1) - (A-1)*cos(w0) - 2*sqrt(A)*alpha

おまけ：バイリニア変換（周波数warpingの補償付き）の変形

                                  1         1 - z^-1
      (normalized)   s  <--  ----------- * ----------
                              tan(w0/2)     1 + z^-1

　それから等式としてこいつらを利用するんだ

                     sin(w0) 
      tan(w0/2) = -------------
                   1 + cos(w0) 
                            1 - cos(w0)
           (tan(w0/2))^2 = -------------
                            1 + cos(w0)

　変形するとこうだ

                 1 + cos(w0)     1 + 2*z^-1 + z^-2
      1    <--  ------------- * -------------------
                 1 + cos(w0)     1 + 2*z^-1 + z^-2

                 1 + cos(w0)     1 - z^-1
      s    <--  ------------- * ----------
                   sin(w0)       1 + z^-1

                  1 + cos(w0)     1         -  z^-2
              =  ------------- * -------------------
                    sin(w0)       1 + 2*z^-1 + z^-2


                 1 + cos(w0)     1 - 2*z^-1 + z^-2
      s^2  <--  ------------- * -------------------
                 1 - cos(w0)     1 + 2*z^-1 + z^-2

　　因数：

                    1 + cos(w0)
                -------------------
                 1 + 2*z^-1 + z^-2

これが全部の分子分母の項で共通なんで、分解できるぜ。すると・・・

                 1 + 2*z^-1 + z^-2
      1    <--  -------------------
                    1 + cos(w0)

                 1         -  z^-2
      s    <--  -------------------
                      sin(w0)


                 1 - 2*z^-1 + z^-2
      s^2  <--  -------------------
                    1 - cos(w0)

さらに、全部の分子分母の項に(sin(w0))^2を掛けると、こうなる。

      1         <--   (1 + 2*z^-1 + z^-2) * (1 - cos(w0))

      s         <--   (1         -  z^-2) * sin(w0)

      s^2       <--   (1 - 2*z^-1 + z^-2) * (1 + cos(w0))

      1 + s^2   <--   2 * (1 - 2*cos(w0)*z^-1 + z^-2)

Biquadの係数は上の式をちょっと簡単にしたら出てくるぜ！

VST の作り方
04.フィルタ

4.1 フィルタの料理法

サイト内検索