Faust : biquad filter用のライブラリを作ってみた

[Faust : 初めてのVSTプラグインを直列にしてみた](https://qiita.com/ttatsf/items/b5fa37babb8ae6f26462)の続き。
FaustLiveのAPFの例に

```Java:APF.dsp
// WARNING: This a "legacy example based on a deprecated library". Check filters.lib
// for more accurate examples of filter functions
```
と書いてあったので ```filters.lib``` を眺めてみたが、似たようなものが見つけられず、とりあえず ```fi.tf2``` を使えっていうことなのかな? ということだけ理解しました。(いや本当は違うのかもしれない...)

係数をどうしていいか分からないので、ググッたり、数式丸パクったりしました。
参考にした記事：

+ [ Cookbook formulae for audio EQ biquad filter coefficients](http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt)
+ [BiQuadフィルタの料理法](https://www.g200kg.com/jp/docs/makingeffects/78743dea3f70c8c2f081b7d5187402ec75e6a6b8.html)

欲しいのはAPFだけだったのですが、どうせならみんな使えるようになると後々幸せになれるのかなと思い、全部入りのを作ってみました。

```Java:RBJ.lib
declare name "RBJ Cookbook Library ";
import("stdfaust.lib");

/////////////////////////////////////////////////////////////////////
// An implementation of biquad filters
// from "Cookbook formulae for audio EQ biquad filter coefficients"
// by Robert Bristow-Johnson  <rbj@audioimagination.com>
// URL : http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
//////////////////////////////////////////////////////////////////////

LPF( F, Q ) 	= rbj( F, 0, Q ).LPF_coef, _ : biquad;
HPF( F, Q ) 	= rbj( F, 0, Q ).HPF_coef, _ : biquad;
BPF( F, Q ) 	= rbj( F, 0, Q ).BPF_coef, _ : biquad;
notch( F, Q ) 	= rbj( F, 0, Q ).notch_coef, _ : biquad;
APF( F, Q ) 	= rbj( F, 0, Q ).APF_coef, _ : biquad;
peakingEQ( F, G, Q ) 	= rbj( F, G, Q ).peakingEQ_coef, _ : biquad;
lowShelf( F, G, Q ) 	= rbj( F, G, Q ).lowShelf_coef, _ : biquad;
highShelf( F, G, Q ) 	= rbj( F, G, Q ).highShelf_coef, _ : biquad;

biquad  = fi.tf2;

rbj( F, G, Q ) = environment{

  LPF_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      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_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      b0 =     ( 1 + cos( w0 ) ) / 2 ;
      b1 = 0 - ( 1 + cos( w0 ) ) ;
      b2 =     ( 1 + cos( w0 ) ) / 2 ;
      a0 =       1 + alpha ;
      a1 =      -2 * cos( w0 ) ;
      a2 =       1 - alpha ;
    } ;

  BPF_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      b0 =   alpha ;
      b1 =   0 ;
      b2 =  -alpha ;
      a0 =   1 + alpha ;
      a1 =  -2 * cos( w0 ) ;
      a2 =   1 - alpha ;
    } ;

  notch_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      b0 =   1 ;
      b1 =  -2 * cos( w0 ) ;
      b2 =   1 ;
      a0 =   1 + alpha ;
      a1 =  -2 * cos( w0 ) ;
      a2 =   1 - alpha ;
    } ;

  APF_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      b0 =   1 - alpha ;
      b1 =  -2 * cos( w0 ) ;
      b2 =   1 + alpha ;
      a0 =   1 + alpha ;
      a1 =  -2 * cos( w0 ) ;
      a2 =   1 - alpha ;
    } ;

  peakingEQ_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      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_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      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_coef = b0/a0, b1/a0, b2/a0, a1/a0, a2/a0
    with{
      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 ;
    } ;

  A = 10 ^ ( G / 40 ) ;
  w0 = 2 * ma.PI * F / ma.SR ;
  alpha = sin( w0 ) / ( 2 * max( 0.001, Q )) ;

};
```
これに合わせて呼び出し側のファイルも変更しました。

```Java:APFDelaySeq.dsp
declare name "APFDelaySeq";
import("stdfaust.lib");
import("RBJ.lib");

process=
  _, _, hslider( "Mid Delay[unit:ms]", 2, 2, 24, 2 )
  : seq( i, 12, bypassIfZero( APFDelay ) )
  : _, _, !
;

bypassIfZero( func, x, y, n ) =
  ( x, y : ba.bypass2( n==0, func ) ), max( 0, n-2 )
;

APFDelay =
  APFSt( 250, 1.5 )
  : APFSt( 700, 4.5 )
  : APFSt( 1250, 6.3 )
  : APFSt( 1600, 15 )
;

APFSt( F, Q ) =
  APF( F, Q ), APF( F, Q )
;
```

ちょっと試してみた感触としては、ちゃんと動いているようです。
おかしいよ、とかあればコメントおねがいします。
## 追記：こっちの方がわかりやすいかも
どのフィルタも同じようなことをしているのでまとめて書こうとしたら、かえって難しい書き方になってしまったのかもです。
上と同じですが、できるだけ簡単な書きかたをしてみました。

```Java:RBJ.lib
declare name "RBJ Cookbook Library ";
import("stdfaust.lib");

/////////////////////////////////////////////////////////////////////
// An implementation of biquad filters
// from "Cookbook formulae for audio EQ biquad filter coefficients"
// by Robert Bristow-Johnson  <rbj@audioimagination.com>
// URL : http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
//////////////////////////////////////////////////////////////////////

LPF( F, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    b0 =  ( 1 - cos( w0 ) ) / 2 ;
    b1 =    1 - cos( w0 ) ;
    b2 =  ( 1 - cos( w0 ) ) / 2 ;
    a0 =    1 + alpha ;
    a1 =   -2 * cos( w0 ) ;
    a2 =    1 - alpha ;

    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

HPF( F, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    b0 =     ( 1 + cos( w0 ) ) / 2 ;
    b1 = 0 - ( 1 + cos( w0 ) ) ;
    b2 =     ( 1 + cos( w0 ) ) / 2 ;
    a0 =       1 + alpha ;
    a1 =      -2 * cos( w0 ) ;
    a2 =       1 - alpha ;

    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

BPF( F, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    b0 =   alpha ;
    b1 =   0 ;
    b2 =  -alpha ;
    a0 =   1 + alpha ;
    a1 =  -2 * cos( w0 ) ;
    a2 =   1 - alpha ;

    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

notch( F, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    b0 =   1 ;
    b1 =  -2 * cos( w0 ) ;
    b2 =   1 ;
    a0 =   1 + alpha ;
    a1 =  -2 * cos( w0 ) ;
    a2 =   1 - alpha ;

    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

APF( F, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    b0 =   1 - alpha ;
    b1 =  -2 * cos( w0 ) ;
    b2 =   1 + alpha ;
    a0 =   1 + alpha ;
    a1 =  -2 * cos( w0 ) ;
    a2 =   1 - alpha ;

    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

peakingEQ( F, G, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    b0 =   1 + alpha * A ;
    b1 =  -2 * cos( w0 ) ;
    b2 =   1 - alpha * A ;
    a0 =   1 + alpha / A ;
    a1 =  -2 * cos( w0 ) ;
    a2 =   1 - alpha / A ;

    A = f_A( G ) ;
    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

lowShelf( F, G, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    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 ;

    A = f_A( G ) ;
    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

highShelf( F, G, Q ) = biquad( b0/a0, b1/a0, b2/a0, a1/a0, a2/a0 )
  with{
    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 ;

    A = f_A( G ) ;
    w0 = f_w0( F ) ;
    alpha = f_alpha( w0, Q ) ;
  } ;

biquad  = fi.tf2;

f_A( G ) = 10 ^ ( G / 40 ) ;
f_w0( F ) = 2 * ma.PI * F / ma.SR ;
f_alpha( w0, Q ) = sin( w0 ) / ( 2 * max( 0.001, Q )) ;
```