What's wrong with this implementation of my IIR filter?

Question

I posted this originally in the electronics stack exchange, as it's for an electronic purpose (digitally controlled boost) but it was pointed out to me that this would be the more appropriate place - so here goes!

As mentioned, I recently started developing a digital regulator for a DCM boost and I've implemented the following discrete transfer function in C:

$$ G_{Z} = \frac{1.39 z^{-1} - 1.332z^{-2} }{2 - 2.002{z^{-1}} + 0.002031z^{-2}} $$

and, in terms of a difference equation:

$$ y[n] = 0.695\cdot y[n-1] - 0.666 \cdot y[n-2] + 1.001\cdot x[n-1] - 0.0010155 \cdot x[n-2] $$

My code is as follows:

// Digi Reg Coeffs
const float a1 = 1.001;
const float a2 = 0.0010155;

const float b1 = 0.6975;
const float b2 = 0.666;

typedef struct {
  float y0; // new reg value
  float y1; // last reg value
  float y2; // last^2 reg value

  float x0; // latest sample
  float x1; // last sample
  float x2; // last^2 sample
} voltage_control;

uint16_t DigitalRegulator(float error) {
  float duty = 0;

  voltage_control.x2 = voltage_control.x1;
  voltage_control.x1 = voltage_control.x0;
  voltage_control.x0 = error; /// this is the new sample

  voltage_control.y2 = voltage_control.y1;
  voltage_control.y1 = voltage_control.y0;

  // calc new y0
  voltage_control.y0 = (a1 * voltage_control.y1) - (a2 * voltage_control.y2) +
                       (b1 * voltage_control.x1) - (b2 * voltage_control.x2);

  // new value for PWM
  duty = (voltage_control.y0 / PWM_GAIN) * PERIOD_VALUE;

  return duty;
}

Is this the correct way to implement a difference equation? I'm asking because the resulting controller seems unstable, and I'd like to know if it's my implementation or something else.

Answer 1

If you take the Z transform on both sides of the difference equation:

$$Y(z)(1-0.695z^{-1}+0.666z^{-2})=X(z)(1.001z^{-1}-0.0010155z^{-2})$$

$$G(z) = \frac{Y(z)}{X(z)}=\frac{1.001z^{-1}-0.0010155z^{-2}}{1-0.695z^{-1}+0.666z^{-2}}=\frac{2.002z^{-1}-0.002031z^{-2}}{2-1.39z^{-1}+1.332z^{-2}}$$

It does not match your transfer function. I think you can work backwards to get the right difference equation now (note the x coefficients go on the numerator and the y's on the denominator). The system you implemented has a pole outside the unit circle so it's unstable.

Answer 2

For your transfer function, the difference equation will be, $$2y[n]-2.002y[n-1]+0.002031y[n-2] = 1.39x[n-1]-1.332x[n-2]$$ $$y[n] = 1.001y[n-1]-0.0010155y[n-2]+0.695x[n-1]-0.666x[n-2]$$

You will have to flip the coefficients.