Biquad フィルタの比較

template<typename Sample> struct ExpSmoother {
  Sample value = 0;
  Sample target = 0;
  Sample kp = 1; // Exponential moving average coefficient.

  void setCutoff(Sample sampleRate, Sample cutoffHz)
  {
    // `double` is used for accuracy.
    double y = double(1)
      - std::cos(double(2) * std::numbers::pi_v<double> * cutoffHz / sampleRate);
    kp = Sample(std::sqrt((y + double(2)) * y) - y);
  }

  void reset(Sample newTarget = 0)
  {
    value = newTarget;
    target = newTarget;
  }

  void push(Sample newTarget = 0) { target = newTarget; }
  inline Sample process() { return value += kp * (target - value); }
};

#include <array>
#include <cmath>
#include <numbers>

template<typename Sample> class RbjBiquad {
  Sample x1 = 0;
  Sample x2 = 0;
  Sample y1 = 0;
  Sample y2 = 0;
  std::array<ExpSmoother<Sample>, 5> co; // Coefficients {b0, b1, b2, -a1, -a2}.

public:
  const std::string name{"RbjBiquad"};

  RbjBiquad()
  {
    for (auto &x : co) x.setCutoff(Sample(1), Sample(0.02));
  }

#define ASSIGN_COEFFICINETS(method)                                                      \
  /* normalizedFreq = cutoffHz / sampleRate. */                                          \
  void method##Lowpass(Sample normalizedFreq, Sample Q)                                  \
  {                                                                                      \
    auto omega = std::numbers::pi_v<Sample> * Sample(2) * normalizedFreq;                \
    auto cos = std::cos(omega);                                                          \
    auto sin = std::sin(omega);                                                          \
    auto alpha = sin / (Sample(2) * Q);                                                  \
    auto a0_inv = Sample(1) / (Sample(1) + alpha);                                       \
    co[0].method(a0_inv *((Sample(1) - cos) / Sample(2))); /* b0 */                      \
    co[1].method(a0_inv *(Sample(1) - cos));               /* b1 */                      \
    co[2].method(a0_inv *((Sample(1) - cos) / Sample(2))); /* b2 */                      \
    co[3].method(-a0_inv *(Sample(-2) * cos));             /* a1 */                      \
    co[4].method(-a0_inv *(Sample(1) - alpha));            /* a2 */                      \
  }

  ASSIGN_COEFFICINETS(reset);
  ASSIGN_COEFFICINETS(push);
#undef ASSIGN_COEFFICINETS

  void reset(Sample normalizedFreq, Sample Q, Sample value = 0)
  {
    x1 = value;
    x2 = value;
    y1 = value;
    y2 = value;

    resetLowpass(normalizedFreq, Q);
  }

  void prepare(Sample normalizedFreq, Sample Q) { pushLowpass(normalizedFreq, Q); }

  Sample process(Sample x0)
  {
    for (auto &x : co) x.process();

    auto y0 = co[0].value * x0 + co[1].value * x1 + co[2].value * x2 + co[3].value * y1
      + co[4].value * y2;

    x2 = x1;
    x1 = x0;
    y2 = y1;
    y1 = y0;

    return y0;
  }
};

template<typename Sample> class TptBiquad {
private:
  ExpSmoother<Sample> g;
  ExpSmoother<Sample> R;
  ExpSmoother<Sample> d;

  Sample s1 = 0;
  Sample s2 = 0;

public:
  const std::string name{"TptBiquad"};

  TptBiquad()
  {
    constexpr Sample cutoff = Sample(0.02);

    g.setCutoff(Sample(1), cutoff);
    R.setCutoff(Sample(1), cutoff);
    d.setCutoff(Sample(1), cutoff);
  }

#define ASSIGN_COEFFICINETS(method)                                                      \
  /* normalizedFreq = cutoffHz / sampleRate. */                                          \
  void method##Coefficient(Sample normalizedFreq, Sample Q)                              \
  {                                                                                      \
    g.method(std::tan(std::numbers::pi_v<Sample> *normalizedFreq));                      \
    R.method(Sample(1) / Q);                                                             \
    d.method(Sample(1) / (Sample(1) + g.target * g.target + g.target * R.target));       \
  }

  ASSIGN_COEFFICINETS(push);
  ASSIGN_COEFFICINETS(reset);
#undef ASSIGN_COEFFICINETS

  void reset(Sample normalizedFreq, Sample Q, Sample value = 0)
  {
    s1 = value;
    s2 = value;
    resetCoefficient(normalizedFreq, Q);
  }

  void prepare(Sample normalizedFreq, Sample Q) { pushCoefficient(normalizedFreq, Q); }

  // Lowpass.
  Sample process(Sample x0)
  {
    g.process();
    d.process();

    auto v1 = (s1 + g.value * (x0 - s2)) * d.value;
    auto v2 = s2 + g.value * v1;
    s1 = Sample(2) * v1 - s1;
    s2 = Sample(2) * v2 - s2;

    return v2;
  }
};