alc/filters/biquad.cpp

   1
   2 #include "config.h"
   3
   4 #include "biquad.h"
   5
   6 #include <algorithm>
   7 #include <cassert>
   8 #include <cmath>
   9
  10 #include "opthelpers.h"
  11
  12
  13 template<typename Real>
  14 void BiquadFilterR<Real>::setParams(BiquadType type, Real gain, Real f0norm, Real rcpQ)
  15 {
  16     // Limit gain to -100dB
  17     assert(gain > 0.00001f);
  18
  19     const Real w0{al::MathDefs<Real>::Tau() * f0norm};
  20     const Real sin_w0{std::sin(w0)};
  21     const Real cos_w0{std::cos(w0)};
  22     const Real alpha{sin_w0/2.0f * rcpQ};
  23
  24     Real sqrtgain_alpha_2;
  25     Real a[3]{ 1.0f, 0.0f, 0.0f };
  26     Real b[3]{ 1.0f, 0.0f, 0.0f };
  27
  28     /* Calculate filter coefficients depending on filter type */
  29     switch(type)
  30     {
  31         case BiquadType::HighShelf:
  32             sqrtgain_alpha_2 = 2.0f * std::sqrt(gain) * alpha;
  33             b[0] =       gain*((gain+1.0f) + (gain-1.0f)*cos_w0 + sqrtgain_alpha_2);
  34             b[1] = -2.0f*gain*((gain-1.0f) + (gain+1.0f)*cos_w0                   );
  35             b[2] =       gain*((gain+1.0f) + (gain-1.0f)*cos_w0 - sqrtgain_alpha_2);
  36             a[0] =             (gain+1.0f) - (gain-1.0f)*cos_w0 + sqrtgain_alpha_2;
  37             a[1] =  2.0f*     ((gain-1.0f) - (gain+1.0f)*cos_w0                   );
  38             a[2] =             (gain+1.0f) - (gain-1.0f)*cos_w0 - sqrtgain_alpha_2;
  39             break;
  40         case BiquadType::LowShelf:
  41             sqrtgain_alpha_2 = 2.0f * std::sqrt(gain) * alpha;
  42             b[0] =       gain*((gain+1.0f) - (gain-1.0f)*cos_w0 + sqrtgain_alpha_2);
  43             b[1] =  2.0f*gain*((gain-1.0f) - (gain+1.0f)*cos_w0                   );
  44             b[2] =       gain*((gain+1.0f) - (gain-1.0f)*cos_w0 - sqrtgain_alpha_2);
  45             a[0] =             (gain+1.0f) + (gain-1.0f)*cos_w0 + sqrtgain_alpha_2;
  46             a[1] = -2.0f*     ((gain-1.0f) + (gain+1.0f)*cos_w0                   );
  47             a[2] =             (gain+1.0f) + (gain-1.0f)*cos_w0 - sqrtgain_alpha_2;
  48             break;
  49         case BiquadType::Peaking:
  50             b[0] =  1.0f + alpha * gain;
  51             b[1] = -2.0f * cos_w0;
  52             b[2] =  1.0f - alpha * gain;
  53             a[0] =  1.0f + alpha / gain;
  54             a[1] = -2.0f * cos_w0;
  55             a[2] =  1.0f - alpha / gain;
  56             break;
  57
  58         case BiquadType::LowPass:
  59             b[0] = (1.0f - cos_w0) / 2.0f;
  60             b[1] =  1.0f - cos_w0;
  61             b[2] = (1.0f - cos_w0) / 2.0f;
  62             a[0] =  1.0f + alpha;
  63             a[1] = -2.0f * cos_w0;
  64             a[2] =  1.0f - alpha;
  65             break;
  66         case BiquadType::HighPass:
  67             b[0] =  (1.0f + cos_w0) / 2.0f;
  68             b[1] = -(1.0f + cos_w0);
  69             b[2] =  (1.0f + cos_w0) / 2.0f;
  70             a[0] =   1.0f + alpha;
  71             a[1] =  -2.0f * cos_w0;
  72             a[2] =   1.0f - alpha;
  73             break;
  74         case BiquadType::BandPass:
  75             b[0] =  alpha;
  76             b[1] =  0.0f;
  77             b[2] = -alpha;
  78             a[0] =  1.0f + alpha;
  79             a[1] = -2.0f * cos_w0;
  80             a[2] =  1.0f - alpha;
  81             break;
  82     }
  83
  84     mA1 = a[1] / a[0];
  85     mA2 = a[2] / a[0];
  86     mB0 = b[0] / a[0];
  87     mB1 = b[1] / a[0];
  88     mB2 = b[2] / a[0];
  89 }
  90
  91 template<typename Real>
  92 void BiquadFilterR<Real>::process(Real *dst, const Real *src, const size_t numsamples)
  93 {
  94     ASSUME(numsamples > 0);
  95
  96     const Real b0{mB0};
  97     const Real b1{mB1};
  98     const Real b2{mB2};
  99     const Real a1{mA1};
 100     const Real a2{mA2};
 101     Real z1{mZ1};
 102     Real z2{mZ2};
 103
 104     /* Processing loop is Transposed Direct Form II. This requires less storage
 105      * compared to Direct Form I (only two delay components, instead of a four-
 106      * sample history; the last two inputs and outputs), and works better for
 107      * floating-point which favors summing similarly-sized values while being
 108      * less bothered by overflow.
 109      *
 110      * See: http://www.earlevel.com/main/2003/02/28/biquads/
 111      */
 112     auto proc_sample = [b0,b1,b2,a1,a2,&z1,&z2](Real input) noexcept -> Real
 113     {
 114         Real output = input*b0 + z1;
 115         z1 = input*b1 - output*a1 + z2;
 116         z2 = input*b2 - output*a2;
 117         return output;
 118     };
 119     std::transform(src, src+numsamples, dst, proc_sample);
 120
 121     mZ1 = z1;
 122     mZ2 = z2;
 123 }
 124
 125 template class BiquadFilterR<float>;
 126 template class BiquadFilterR<double>;