Filter.h

Go to the documentation of this file.

#pragma once
#include <math.h>
 
#include "AudioToolsConfig.h"
#include "AudioTools/CoreAudio/AudioLogger.h"
#ifdef USE_TYPETRAITS
#include <type_traits>
#endif
 
 
namespace audio_tools {
 
template <typename T = float>
class Filter {
 public:
  // construct without coefs
  Filter() = default;
  virtual ~Filter() = default;
  Filter(Filter const&) = delete;
  Filter& operator=(Filter const&) = delete;
  virtual T process(T in) = 0;
  virtual void reset() {}
};
 
template <typename T = float>
class NoFilter : public Filter<T> {
 public:
  // construct without coefs
  NoFilter() = default;
  T process(T in) override { return in; }
};
 
template <typename T = float>
class FIR : public Filter<T> {
 public:
  template <size_t B>
  FIR(const T (&b)[B], const T factor = 1.0)
      : lenB(B), factor(factor), needs_divide(factor != T(1)) {
    setValues(b);
  }
 
  template <size_t B>
  void setValues(const T (&b)[B]) {
    x.resize(lenB);
    coeff_b.resize(2 * lenB - 1);
    for (uint16_t i = 0; i < 2 * lenB - 1; i++) {
      coeff_b[i] = b[(2 * lenB - 1 - i) % lenB];
    }
  }
 
  void reset() override {
    i_b = 0;
    for (uint16_t i = 0; i < lenB; i++) x[i] = 0;
  }
 
  T process(T value) override {
    x[i_b] = value;
    T b_terms = 0;
    T* b_shift = &coeff_b[lenB - i_b - 1];
    for (uint16_t i = 0; i < lenB; i++) {
      b_terms += b_shift[i] * x[i];
    }
    i_b++;
    if (i_b == lenB) i_b = 0;
 
#ifdef USE_TYPETRAITS
    if (!(std::is_same<T, float>::value || std::is_same<T, double>::value)) {
      b_terms = b_terms / factor;
    }
#else
    // needs_divide is computed once in the constructor (not per sample) --
    // for a bounded fixed-point T like q1_14_t, comparing against a float
    // literal here on every process() call would otherwise force a
    // per-sample float construction/comparison.
    if (needs_divide) {
      b_terms = b_terms / factor;
    }
#endif
    return b_terms;
  }
 
 private:
  const uint16_t lenB;
  uint16_t i_b = 0;
  Vector<T> x;
  Vector<T> coeff_b;
  T factor;
  bool needs_divide;
};
 
template <typename T = float>
class IIR : public Filter<T> {
 public:
  template <size_t B, size_t A>
  IIR(const T (&b)[B], const T (&_a)[A], T factor = 1.0)
      : factor(factor), needs_divide(factor != T(1)), lenB(B), lenA(A - 1) {
    x.resize(lenB);
    y.resize(lenA);
    coeff_b.resize(2 * lenB - 1);
    coeff_a.resize(2 * lenA - 1);
    T a0 = _a[0];
    const T* a = &_a[1];
    for (uint16_t i = 0; i < 2 * lenB - 1; i++) {
      coeff_b[i] = b[(2 * lenB - 1 - i) % lenB] / a0;
    }
    for (uint16_t i = 0; i < 2 * lenA - 1; i++) {
      coeff_a[i] = a[(2 * lenA - 2 - i) % lenA] / a0;
    }
  }
 
  void reset() override {
    i_b = 0;
    i_a = 0;
    for (uint16_t i = 0; i < lenB; i++) x[i] = 0;
    for (uint16_t i = 0; i < lenA; i++) y[i] = 0;
  }
 
  T process(T value) override {
    x[i_b] = value;
    T b_terms = 0;
    T* b_shift = &coeff_b[lenB - i_b - 1];
 
    T a_terms = 0;
    T* a_shift = &coeff_a[lenA - i_a - 1];
 
    for (uint16_t i = 0; i < lenB; i++) {
      b_terms += x[i] * b_shift[i];
    }
    for (uint16_t i = 0; i < lenA; i++) {
      a_terms += y[i] * a_shift[i];
    }
 
    T filtered = b_terms - a_terms;
    y[i_a] = filtered;
    i_b++;
    if (i_b == lenB) i_b = 0;
    i_a++;
    if (i_a == lenA) i_a = 0;
 
#ifdef USE_TYPETRAITS
    if (!(std::is_same<T, float>::value || std::is_same<T, double>::value)) {
      filtered = filtered / factor;
    }
#else
    // needs_divide is computed once in the constructor (not per sample) --
    // see FIR::process() for why.
    if (needs_divide) {
      filtered = filtered / factor;
    }
#endif
    return filtered;
  }
 
 private:
  T factor;
  bool needs_divide;
  const uint16_t lenB, lenA;
  uint16_t i_b = 0, i_a = 0;
  Vector<T> x;
  Vector<T> y;
  Vector<T> coeff_b;
  Vector<T> coeff_a;
};
 
template <typename T = float>
class BiQuadDF1 : public Filter<T> {
 public:
  BiQuadDF1(const T (&b)[3], const T (&a)[3])
      : b_0(b[0] / a[0]),
        b_1(b[1] / a[0]),
        b_2(b[2] / a[0]),
        a_1(a[1] / a[0]),
        a_2(a[2] / a[0]) {}
  BiQuadDF1(const T (&b)[3], const T (&a)[2])
      : b_0(b[0]), b_1(b[1]), b_2(b[2]), a_1(a[0]), a_2(a[1]) {}
  BiQuadDF1(const T (&b)[3], const T (&a)[2], T gain)
      : b_0(gain * b[0]),
        b_1(gain * b[1]),
        b_2(gain * b[2]),
        a_1(a[0]),
        a_2(a[1]) {}
  BiQuadDF1(const T (&b)[3], const T (&a)[3], T gain)
      : b_0(gain * b[0] / a[0]),
        b_1(gain * b[1] / a[0]),
        b_2(gain * b[2] / a[0]),
        a_1(a[1] / a[0]),
        a_2(a[2] / a[0]) {}
 
  void reset() override { x_0 = x_1 = y_1 = y_2 = 0; }
 
  T process(T value) override {
    T x_2 = x_1;
    x_1 = x_0;
    x_0 = value;
    T b_terms = x_0 * b_0 + x_1 * b_1 + x_2 * b_2;
    T a_terms = y_1 * a_1 + y_2 * a_2;
    y_2 = y_1;
    y_1 = b_terms - a_terms;
    return y_1;
  }
 
 protected:
  T b_0 = 0;
  T b_1 = 0;
  T b_2 = 0;
  T a_1 = 0;
  T a_2 = 0;
 
  // allow constructor w/o parameter in subclasses
  BiQuadDF1() = default;
 
  T x_0 = 0;
  T x_1 = 0;
  T y_1 = 0;
  T y_2 = 0;
};
 
template <typename T = float>
class BiQuadDF2 : public Filter<T> {
 public:
  BiQuadDF2(const T (&b)[3], const T (&a)[3])
      : b_0(b[0] / a[0]),
        b_1(b[1] / a[0]),
        b_2(b[2] / a[0]),
        a_1(a[1] / a[0]),
        a_2(a[2] / a[0]) {}
  BiQuadDF2(const T (&b)[3], const T (&a)[2])
      : b_0(b[0]), b_1(b[1]), b_2(b[2]), a_1(a[0]), a_2(a[1]) {}
  BiQuadDF2(const T (&b)[3], const T (&a)[2], T gain)
      : b_0(gain * b[0]),
        b_1(gain * b[1]),
        b_2(gain * b[2]),
        a_1(a[0]),
        a_2(a[1]) {}
  BiQuadDF2(const T (&b)[3], const T (&a)[3], T gain)
      : b_0(gain * b[0] / a[0]),
        b_1(gain * b[1] / a[0]),
        b_2(gain * b[2] / a[0]),
        a_1(a[1] / a[0]),
        a_2(a[2] / a[0]) {}
 
  void reset() override { w_0 = w_1 = 0; }
 
  T process(T value) override {
    T w_2 = w_1;
    w_1 = w_0;
    w_0 = value - a_1 * w_1 - a_2 * w_2;
    T y = b_0 * w_0 + b_1 * w_1 + b_2 * w_2;
    return y;
  }
 
 protected:
  T b_0 = 0;
  T b_1 = 0;
  T b_2 = 0;
  T a_1 = 0;
  T a_2 = 0;
 
  // allow constructor w/o parameter in subclasses
  BiQuadDF2() = default;
 
  T w_0 = 0;
  T w_1 = 0;
};
 
struct BiQuadCoeffs {
  float b_0, b_1, b_2, a_1, a_2;
};
 
template <typename T>
void checkCoeffRange(const BiQuadCoeffs &c) {
  auto check = [](const char *name, float value) {
    T t = value;
    float roundtrip = (float)t;
    // allow normal quantization noise (e.g. ~1/16384 for q1_14_t); anything
    // larger indicates the value was clamped/saturated, not just rounded.
    if (fabs(roundtrip - value) > fabs(value) * 0.01f + 1e-4f) {
      LOGE(
          "Filter coefficient %s=%f does not fit in the range of T (stored "
          "as %f) - the filter will be inaccurate",
          name, value, roundtrip);
    }
  };
  check("b_0", c.b_0);
  check("b_1", c.b_1);
  check("b_2", c.b_2);
  check("a_1", c.a_1);
  check("a_2", c.a_2);
}
 
inline BiQuadCoeffs calculateLowPassCoeffs(float frequency, float sampleRate,
                                            float q) {
  float w0 = frequency * (2.0f * PI / sampleRate);
  float sinW0 = sin(w0);
  float alpha = sinW0 / (q * 2.0f);
  float cosW0 = cos(w0);
  float scale = 1.0f / (1.0f + alpha);
  BiQuadCoeffs c;
  c.b_0 = ((1.0f - cosW0) / 2.0f) * scale;
  c.b_1 = (1.0f - cosW0) * scale;
  c.b_2 = c.b_0;
  c.a_1 = (-2.0f * cosW0) * scale;
  c.a_2 = (1.0f - alpha) * scale;
  return c;
}
 
template <typename T = float>
class LowPassFilter : public BiQuadDF2<T> {
 public:
  LowPassFilter() = default;
  LowPassFilter(float frequency, float sampleRate, float q = 0.7071f)
      : BiQuadDF2<T>() {
    begin(frequency, sampleRate, q);
  }
  void begin(float frequency, float sampleRate, float q = 0.7071f) {
    BiQuadCoeffs c = calculateLowPassCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF2<T>::b_0 = c.b_0;
    BiQuadDF2<T>::b_1 = c.b_1;
    BiQuadDF2<T>::b_2 = c.b_2;
    BiQuadDF2<T>::a_1 = c.a_1;
    BiQuadDF2<T>::a_2 = c.a_2;
  }
};
 
template <typename T = float>
class LowPassFilterDF1 : public BiQuadDF1<T> {
 public:
  LowPassFilterDF1() = default;
  LowPassFilterDF1(float frequency, float sampleRate, float q = 0.7071f)
      : BiQuadDF1<T>() {
    begin(frequency, sampleRate, q);
  }
  void begin(float frequency, float sampleRate, float q = 0.7071f) {
    BiQuadCoeffs c = calculateLowPassCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF1<T>::b_0 = c.b_0;
    BiQuadDF1<T>::b_1 = c.b_1;
    BiQuadDF1<T>::b_2 = c.b_2;
    BiQuadDF1<T>::a_1 = c.a_1;
    BiQuadDF1<T>::a_2 = c.a_2;
  }
};
 
inline BiQuadCoeffs calculateHighPassCoeffs(float frequency, float sampleRate,
                                             float q) {
  float w0 = frequency * (2.0f * PI / sampleRate);
  float sinW0 = sin(w0);
  float alpha = sinW0 / (q * 2.0f);
  float cosW0 = cos(w0);
  float scale = 1.0f / (1.0f + alpha);
  BiQuadCoeffs c;
  c.b_0 = ((1.0f + cosW0) / 2.0f) * scale;
  c.b_1 = -(1.0f + cosW0) * scale;
  c.b_2 = c.b_0;
  c.a_1 = (-2.0f * cosW0) * scale;
  c.a_2 = (1.0f - alpha) * scale;
  return c;
}
 
template <typename T = float>
class HighPassFilter : public BiQuadDF2<T> {
 public:
  HighPassFilter() = default;
  HighPassFilter(float frequency, float sampleRate, float q = 0.7071f)
      : BiQuadDF2<T>() {
    begin(frequency, sampleRate, q);
  }
  void begin(float frequency, float sampleRate, float q = 0.7071f) {
    BiQuadCoeffs c = calculateHighPassCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF2<T>::b_0 = c.b_0;
    BiQuadDF2<T>::b_1 = c.b_1;
    BiQuadDF2<T>::b_2 = c.b_2;
    BiQuadDF2<T>::a_1 = c.a_1;
    BiQuadDF2<T>::a_2 = c.a_2;
  }
};
 
template <typename T = float>
class HighPassFilterDF1 : public BiQuadDF1<T> {
 public:
  HighPassFilterDF1() = default;
  HighPassFilterDF1(float frequency, float sampleRate, float q = 0.7071f)
      : BiQuadDF1<T>() {
    begin(frequency, sampleRate, q);
  }
  void begin(float frequency, float sampleRate, float q = 0.7071f) {
    BiQuadCoeffs c = calculateHighPassCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF1<T>::b_0 = c.b_0;
    BiQuadDF1<T>::b_1 = c.b_1;
    BiQuadDF1<T>::b_2 = c.b_2;
    BiQuadDF1<T>::a_1 = c.a_1;
    BiQuadDF1<T>::a_2 = c.a_2;
  }
};
 
inline BiQuadCoeffs calculateBandPassCoeffs(float frequency, float sampleRate,
                                             float q) {
  float w0 = frequency * (2.0f * PI / sampleRate);
  float sinW0 = sin(w0);
  float alpha = sinW0 / (q * 2.0f);
  float cosW0 = cos(w0);
  float scale = 1.0f / (1.0f + alpha);
  BiQuadCoeffs c;
  c.b_0 = alpha * scale;
  c.b_1 = 0;
  c.b_2 = (-alpha) * scale;
  c.a_1 = (-2.0f * cosW0) * scale;
  c.a_2 = (1.0f - alpha) * scale;
  return c;
}
 
template <typename T = float>
class BandPassFilter : public BiQuadDF2<T> {
 public:
  BandPassFilter() = default;
  BandPassFilter(float frequency, float sampleRate, float q = 1.0)
      : BiQuadDF2<T>() {
    begin(frequency, sampleRate, q);
  }
  void begin(float frequency, float sampleRate, float q = 1.0) {
    BiQuadCoeffs c = calculateBandPassCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF2<T>::b_0 = c.b_0;
    BiQuadDF2<T>::b_1 = c.b_1;
    BiQuadDF2<T>::b_2 = c.b_2;
    BiQuadDF2<T>::a_1 = c.a_1;
    BiQuadDF2<T>::a_2 = c.a_2;
  }
};
 
template <typename T = float>
class BandPassFilterDF1 : public BiQuadDF1<T> {
 public:
  BandPassFilterDF1() = default;
  BandPassFilterDF1(float frequency, float sampleRate, float q = 1.0)
      : BiQuadDF1<T>() {
    begin(frequency, sampleRate, q);
  }
  void begin(float frequency, float sampleRate, float q = 1.0) {
    BiQuadCoeffs c = calculateBandPassCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF1<T>::b_0 = c.b_0;
    BiQuadDF1<T>::b_1 = c.b_1;
    BiQuadDF1<T>::b_2 = c.b_2;
    BiQuadDF1<T>::a_1 = c.a_1;
    BiQuadDF1<T>::a_2 = c.a_2;
  }
};
 
 
inline BiQuadCoeffs calculateNotchCoeffs(float frequency, float sampleRate,
                                          float q) {
  float w0 = frequency * (2.0f * PI / sampleRate);
  float sinW0 = sin(w0);
  float alpha = sinW0 / (q * 2.0f);
  float cosW0 = cos(w0);
  float scale = 1.0f / (1.0f + alpha);
  BiQuadCoeffs c;
  c.b_0 = scale;
  c.b_1 = (-2.0f * cosW0) * scale;
  c.b_2 = c.b_0;
  c.a_1 = (-2.0f * cosW0) * scale;
  c.a_2 = (1.0f - alpha) * scale;
  return c;
}
 
template <typename T = float>
class NotchFilter : public BiQuadDF2<T> {
 public:
  NotchFilter() = default;
  NotchFilter(float frequency, float sampleRate, float q = 1.0)
      : BiQuadDF2<T>() {
    begin(frequency, sampleRate, q);
  }
 
  void begin(float frequency, float sampleRate, float q = 1.0) {
    BiQuadCoeffs c = calculateNotchCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF2<T>::b_0 = c.b_0;
    BiQuadDF2<T>::b_1 = c.b_1;
    BiQuadDF2<T>::b_2 = c.b_2;
    BiQuadDF2<T>::a_1 = c.a_1;
    BiQuadDF2<T>::a_2 = c.a_2;
  }
};
 
template <typename T = float>
class NotchFilterDF1 : public BiQuadDF1<T> {
 public:
  NotchFilterDF1() = default;
  NotchFilterDF1(float frequency, float sampleRate, float q = 1.0)
      : BiQuadDF1<T>() {
    begin(frequency, sampleRate, q);
  }
 
  void begin(float frequency, float sampleRate, float q = 1.0) {
    BiQuadCoeffs c = calculateNotchCoeffs(frequency, sampleRate, q);
    checkCoeffRange<T>(c);
    BiQuadDF1<T>::b_0 = c.b_0;
    BiQuadDF1<T>::b_1 = c.b_1;
    BiQuadDF1<T>::b_2 = c.b_2;
    BiQuadDF1<T>::a_1 = c.a_1;
    BiQuadDF1<T>::a_2 = c.a_2;
  }
};
 
inline BiQuadCoeffs calculateLowShelfCoeffs(float frequency, float sampleRate,
                                             float gain, float slope) {
  float a = pow(10.0f, gain / 40.0f);
  float w0 = frequency * (2.0f * PI / sampleRate);
  float sinW0 = sin(w0);
  // float alpha = (sinW0 * sqrt((a+1/a)*(1/slope-1)+2) ) / 2.0;
  float cosW0 = cos(w0);
  // generate three helper-values (intermediate results):
  float sinsq =
      sinW0 * sqrt((pow(a, 2.0f) + 1.0f) * (1.0f / slope - 1.0f) + 2.0f * a);
  float aMinus = (a - 1.0f) * cosW0;
  float aPlus = (a + 1.0f) * cosW0;
  float scale = 1.0f / ((a + 1.0f) + aMinus + sinsq);
  BiQuadCoeffs c;
  c.b_0 = a * ((a + 1.0f) - aMinus + sinsq) * scale;
  c.b_1 = 2.0f * a * ((a - 1.0f) - aPlus) * scale;
  c.b_2 = a * ((a + 1.0f) - aMinus - sinsq) * scale;
  c.a_1 = -2.0f * ((a - 1.0f) + aPlus) * scale;
  c.a_2 = ((a + 1.0f) + aMinus - sinsq) * scale;
  return c;
}
 
template <typename T = float>
class LowShelfFilter : public BiQuadDF2<T> {
 public:
  LowShelfFilter() = default;
  LowShelfFilter(float frequency, float sampleRate, float gain,
                 float slope = 1.0f)
      : BiQuadDF2<T>() {
    begin(frequency, sampleRate, gain, slope);
  }
 
  void begin(float frequency, float sampleRate, float gain,
             float slope = 1.0f) {
    BiQuadCoeffs c = calculateLowShelfCoeffs(frequency, sampleRate, gain, slope);
    checkCoeffRange<T>(c);
    BiQuadDF2<T>::b_0 = c.b_0;
    BiQuadDF2<T>::b_1 = c.b_1;
    BiQuadDF2<T>::b_2 = c.b_2;
    BiQuadDF2<T>::a_1 = c.a_1;
    BiQuadDF2<T>::a_2 = c.a_2;
  }
};
 
template <typename T = float>
class LowShelfFilterDF1 : public BiQuadDF1<T> {
 public:
  LowShelfFilterDF1() = default;
  LowShelfFilterDF1(float frequency, float sampleRate, float gain,
                     float slope = 1.0f)
      : BiQuadDF1<T>() {
    begin(frequency, sampleRate, gain, slope);
  }
 
  void begin(float frequency, float sampleRate, float gain,
             float slope = 1.0f) {
    BiQuadCoeffs c = calculateLowShelfCoeffs(frequency, sampleRate, gain, slope);
    checkCoeffRange<T>(c);
    BiQuadDF1<T>::b_0 = c.b_0;
    BiQuadDF1<T>::b_1 = c.b_1;
    BiQuadDF1<T>::b_2 = c.b_2;
    BiQuadDF1<T>::a_1 = c.a_1;
    BiQuadDF1<T>::a_2 = c.a_2;
  }
};
 
inline BiQuadCoeffs calculateHighShelfCoeffs(float frequency,
                                              float sampleRate, float gain,
                                              float slope) {
  float a = pow(10.0f, gain / 40.0f);
  float w0 = frequency * (2.0f * PI / sampleRate);
  float sinW0 = sin(w0);
  // float alpha = (sinW0 * sqrt((a+1/a)*(1/slope-1)+2) ) / 2.0;
  float cosW0 = cos(w0);
  // generate three helper-values (intermediate results):
  float sinsq =
      sinW0 * sqrt((pow(a, 2.0f) + 1.0f) * (1.0f / slope - 1.0f) + 2.0f * a);
  float aMinus = (a - 1.0f) * cosW0;
  float aPlus = (a + 1.0f) * cosW0;
  float scale = 1.0f / ((a + 1.0f) - aMinus + sinsq);
  BiQuadCoeffs c;
  c.b_0 = a * ((a + 1.0f) + aMinus + sinsq) * scale;
  c.b_1 = -2.0f * a * ((a - 1.0f) + aPlus) * scale;
  c.b_2 = a * ((a + 1.0f) + aMinus - sinsq) * scale;
  c.a_1 = 2.0f * ((a - 1.0f) - aPlus) * scale;
  c.a_2 = ((a + 1.0f) - aMinus - sinsq) * scale;
  return c;
}
 
template <typename T = float>
class HighShelfFilter : public BiQuadDF2<T> {
 public:
  HighShelfFilter() = default;
  HighShelfFilter(float frequency, float sampleRate, float gain,
                  float slope = 1.0f)
      : BiQuadDF2<T>() {
    begin(frequency, sampleRate, gain, slope);
  }
  void begin(float frequency, float sampleRate, float gain,
             float slope = 1.0f) {
    BiQuadCoeffs c = calculateHighShelfCoeffs(frequency, sampleRate, gain, slope);
    checkCoeffRange<T>(c);
    BiQuadDF2<T>::b_0 = c.b_0;
    BiQuadDF2<T>::b_1 = c.b_1;
    BiQuadDF2<T>::b_2 = c.b_2;
    BiQuadDF2<T>::a_1 = c.a_1;
    BiQuadDF2<T>::a_2 = c.a_2;
  }
};
 
template <typename T = float>
class HighShelfFilterDF1 : public BiQuadDF1<T> {
 public:
  HighShelfFilterDF1() = default;
  HighShelfFilterDF1(float frequency, float sampleRate, float gain,
                      float slope = 1.0f)
      : BiQuadDF1<T>() {
    begin(frequency, sampleRate, gain, slope);
  }
  void begin(float frequency, float sampleRate, float gain,
             float slope = 1.0f) {
    BiQuadCoeffs c = calculateHighShelfCoeffs(frequency, sampleRate, gain, slope);
    checkCoeffRange<T>(c);
    BiQuadDF1<T>::b_0 = c.b_0;
    BiQuadDF1<T>::b_1 = c.b_1;
    BiQuadDF1<T>::b_2 = c.b_2;
    BiQuadDF1<T>::a_1 = c.a_1;
    BiQuadDF1<T>::a_2 = c.a_2;
  }
};
 
template <typename T, size_t N>
class SOSFilter : public Filter<T> {
 public:
  SOSFilter(const T (&b)[N][3], const T (&a)[N][3], const T (&gain)[N]) {
    for (size_t i = 0; i < N; i++)
      filters[i] = new BiQuadDF2<T>(b[i], a[i], gain[i]);
  }
  SOSFilter(const T (&sos)[N][6], const T (&gain)[N]) {
    for (size_t i = 0; i < N; i++) {
      T b[3];
      T a[3];
      copy(b, &sos[i][0]);
      copy(a, &sos[i][3]);
      filters[i] = new BiQuadDF2<T>(b, a, gain[i]);
    }
  }
  SOSFilter(const T (&b)[N][3], const T (&a)[N][2], const T (&gain)[N]) {
    for (size_t i = 0; i < N; i++)
      filters[i] = new BiQuadDF2<T>(b[i], a[i], gain[i]);
  }
  SOSFilter(const T (&b)[N][3], const T (&a)[N][2]) {
    for (size_t i = 0; i < N; i++) filters[i] = new BiQuadDF2<T>(b[i], a[i]);
  }
  SOSFilter(const T (&b)[N][3], const T (&a)[N][3]) {
    for (size_t i = 0; i < N; i++) filters[i] = new BiQuadDF2<T>(b[i], a[i]);
  }
  SOSFilter(SOSFilter const&) = delete;
  SOSFilter& operator=(SOSFilter const&) = delete;
  ~SOSFilter() {
    for (size_t i = 0; i < N; i++) delete filters[i];
  }
  void reset() override {
    for (Filter<T>*& filter : filters) filter->reset();
  }
  T process(T value) override {
    for (Filter<T>*& filter : filters) value = filter->process(value);
    return value;
  }
 
 private:
  Filter<T>* filters[N];
  template <size_t M>
  void copy(T (&dest)[M], const T* src) {
    for (size_t i = 0; i < M; i++) dest[i] = src[i];
  }
};
 
template <typename T, size_t N>
class FilterChain : public Filter<T> {
 public:
  FilterChain(Filter<T>* (&&filters)[N]) {
    for (size_t i = 0; i < N; i++) {
      this->filters[i] = filters[i];
    }
  }
 
  void reset() override {
    for (Filter<T>*& filter : filters) {
      if (filter != nullptr) filter->reset();
    }
  }
 
  T process(T value) override {
    for (Filter<T>*& filter : filters) {
      if (filter != nullptr) {
        value = filter->process(value);
      }
    }
    return value;
  }
 
 private:
  Filter<T>* filters[N] = {0};
};
 
}  // namespace audio_tools