// updates the coefficients with the "ZDF" step invariant method.
void update(float _Q, float F) {
  D = _Q;
  // _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
  // pre iterations values (equiv to Chamberlin filter)
  float _a = F * F;
  float tmp = 1 - _a - D * F;
  float _b = F * tmp + F;
  float _c = tmp * tmp - _a;
  // _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
  // 7 iterations to obtain the new a b and c (equiv x32 oversampling)
  // Two sets of vars (_a _b _c and na nb nc) are used to bounce values
  // between iterations to optimize copies.
  float b2;
  b2 = _b * _b;
  na = b2 + _a * (2 - _a);
  nb = _b * (1 + _c - _a);
  nc = _c * _c - b2;
  b2 = nb * nb;
  _a = b2 + na * (2 - na);
  _b = nb * (1 + nc - na);
  _c = nc * nc - b2;
  b2 = _b * _b;
  na = b2 + _a * (2 - _a);
  nb = _b * (1 + _c - _a);
  nc = _c * _c - b2;
  b2 = nb * nb;
  _a = b2 + na * (2 - na);
  _b = nb * (1 + nc - na);
  _c = nc * nc - b2;
  b2 = _b * _b;
  na = b2 + _a * (2 - _a);
  nb = _b * (1 + _c - _a);
  nc = _c * _c - b2;
  b2 = nb * nb;
  _a = b2 + na * (2 - na);
  _b = nb * (1 + nc - na);
  _c = nc * nc - b2;
  b2 = _b * _b;
  na = b2 + _a * (2 - _a);
  nb = _b * (1 + _c - _a);
  nc = _c * _c - b2;
  // _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
  // deltas for sample rate interpolation
  da = (na - a) * (1 / 32.0f);
  db = (nb - b) * (1 / 32.0f);
  dc = (nc - c) * (1 / 32.0f);
}

float a, b, c;    // current values (srate interpolated)
float na, nb, nc; // new values
float da, db, dc; // deltas for srate interpolation
float D;          // used for hp
float x;          // input
float lp, bp;     // state variables/outputs
float TRF_coef;

D = x = lp = bp = 0;
update(1.0f, 0.1f / 32);
a = na;
b = nb;
c = nc;
da = db = dc = 0;
TRF_coef = (820.0f / (1 << 27)) * 2 * 3.1415926535f / (128 * 48000);

float q = __USAT(param_Q + inlet_reso, 27) * (1.0f / (1 << 27));
q = 0.25f + q * (1 + q * q * (18.75f + q * 60)); // 0.25 -> 80

int32_t alpha;
MTOFEXTENDED(param_pitch + inlet_pitch - (27 << 20), alpha);
update(1 / (2 * q), alpha *TRF_coef);

// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
a += da;
b += db;
c += dc;
x = (float)inlet_x1;
float x_lp = x - lp;
outlet_lp1 = (int32_t)(lp += a * x_lp + b * bp);
outlet_hp1 = (int32_t)(x - D * bp - lp);
outlet_bp1 = (int32_t)(bp = b * x_lp + c * bp);
a += da;
b += db;
c += dc;
// _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
x = (float)inlet_x0;
x_lp = x - lp;
outlet_lp0 = (int32_t)(lp += a * x_lp + b * bp);
outlet_hp0 = (int32_t)(x - D * bp - lp);
outlet_bp0 = (int32_t)(bp = b * x_lp + c * bp);