Digital filtering on a microcontroller using integer math

Question

First of all, please be patient with this post, I don't know how to explain my problem in a shorter way. Thank you in advance! :)

I have a question regarding recursive difference equations implemented using integer variables.

Let us suppose that I want to filter an input signal $u(t)$ with a cut-off frequency $f_f=\frac{1}{2\pi T_f}$. In the continuous time domain, the filter transfer function is:

$$G_\mathcal{L}(s)=\frac{1}{1+sT_f}.$$

The equivalent transfer function in the $\mathcal{Z}$-domain that has the same response to a step function $u(t)=\mu(t)$ (in other words, $G_\mathcal{L}(s)$ is discretized using the ZOH discretization method) is defined as follows:

$$G_\mathcal{Z} = \frac{(1-a)z^{-1}}{1-az^{-1}}, \quad a=e^{-\frac{T_s}{T_f}}, \quad 0<a<1,$$

where $T_s$ is the system sample time. The recursive difference equation is:

$$y_n = u_{n-1} + \frac{K}{M} \bigl( y_{n-1} - u_{n-1} \bigr), \quad a=\frac{K}{M}, \quad n,K,M\in\mathbb{Z}, $$

where $n$ is the time index, $y_n \equiv y(n)$ and $u_n \equiv u(n)$. Note that I usually choose $M$ to be of the form $M=2^b$, where $b \in \mathbb{N}$, so that I don't need to use a divide instruction which requires a lot of instruction cycles on a microcontroller, i.e., I can use a bit shift instruction instead, as follows:

y = u + ( K * ( y - u ) ) >> b );

Please note that in this "integer division" operation using the bit shift (>>b), the result is always rounded towards the first lower number, i.e., 4/3=1 and (-4)/3=-2.

Example 1. Let us calculate an evolution of $y$ for input $u(n)=10000\mu(n)$, initial condition $y(0)=0$, and filter parameters $K=311$ and $M=512$.

This difference equation would evolve on a microcontroller as follows:

k     y(k)    u(k)    y(k+1)
0     0       10000   3925
1     3925    10000   6309
2     6309    10000   7758
...   ...     ...     ...
18    9997    10000   9998
19    9998    10000   9998
20    9998    10000   9998

As you can see, the digital filter output never reaches the input value, i.e. the error is $e_\infty = 2$. How can I formalize this error, i.e., is it possible to calculate the $e_\infty$ from $K$ and $M$ only? I've run several simulations of the system, it seems that the error $e_\infty$ doesn't depend on the input signal. Please find attached a diagram of $e_\infty$ vs. $K$ for a fixed $M=2^9$.

The interesting thing here is if we change the initial condition to be greater than $u$, for example, $y(0)=20000$, then the digital filter manages to reach the final value for all $0<K<M$.

I managed to find a solution to this problem, as follows:

z = K * ( y - u )
if (z < 0) {
    Z = Z + c;
}
y = u + (Z >> b);

where c=(2^b)-1. In other words, I round toward 0 for negative numbers.

However, I'm not comfortable using this solution, as I don't fully understand what happens here. Is there any smarter solution than the one with rounding toward 0? Can you suggest some tutorials (articles) on how to implement algorithms using integer math only. I don't want to reinvent hot water here :)

If you managed to reach until the end of this post, once again, I'm sorry for the large post! :)

Answer 1

I've managed to find an answer to my question. The steady-state error $e_\infty$ can be estimated as:

$$e_\infty = \left\lfloor \frac{M-1}{M-K} \right\rfloor, \quad M,K\in\mathbb{Z}.$$

Here is a proof.

Let us recall the difference equation:

$$y_n = u_{n-1} + \frac{K}{M}\left( y_{n-1}-u_{n-1} \right) , \quad u,y,n\in\mathbb{Z},$$

where $M=2^b$, $b\in\mathbb{N}$, and $-M<K<M$ for a system to be stable.

If we write $e_n = y_n - u_n$, then this equation can be written as:

$$e_n = \frac{K}{M}e_{n-1}, \quad e\in\mathbb{Z},$$

since $u_{k-1} = u_k$. Assuming that the system is stable, this equation can be rewritten for a steady state as follows:

$$e_\infty = \left\lfloor \frac{K}{M}e_\infty \right\rfloor, \quad e_\infty\in\mathbb{Z}.$$

Now let us assume that for the decimal number on the RHS of the previous equation we use a floor function, i.e., $\lfloor1.3\rfloor=1$ and $\lfloor-1.3\rfloor=-2$.

If we assume for example that $e_\infty=-2$, then the RHS of the previous equation (assuming that we use a floor function) needs to be:

$$\frac{K}{M} (-2) < (-1) .$$

Let us solve this inequality:

$$\frac{K}{M} \geq \frac{1}{2} \rightarrow K \geq M \cdot \frac{1}{2}$$

For $M=512$ and $e_\infty=-2$, the parameter $K$ is at least $K=256$. Simulations of the system have shown that for $K=256$ the steady-state error is $e_\infty=-1$, and here we predicted it to be $e_\infty=-2$. I'm still not sure what is wrong in this reasoning, but this can be fixed by replacing $\geq$ with $>$ in the previous ineqaulity. The previous equation can be generalized as follows:

$$K > M \cdot \frac{e_\infty + 1}{e_\infty} .$$

Now let us see how this inequality works for $M=512$ and different $e_\infty$:

if e = -1 then K > 0
if e = -2 then K > 256
if e = -3 then K > 341.33
if e = -4 then K > 384
if e = -5 then K > 409.6

This table has been confirmed by simulating the system for different $K$.

Now let us see how this table looks like for positive errors:

if e = +1 then K > 1024    (not a feasible solution)
if e = +2 then K > 768     (not a feasible solution)
if e = +3 then K > 682.66  (not a feasible solution)

In other words, for any $K < M$, there is no positive steady-state error. This has been also confirmed by simulations. Following the same explanation, if we use a ceil function instead of a floor function, we get the same problem for the positive steady-state errors.

Conclusion: a floor function should be used for positive decimal numbers, and a ceil function should be used for negative decimal numbers. In that case there will be no error in the steady-state condition.

Z = K * ( y - u );
if ( Z < 0 ) {
    Z = Z + ( (1<<b) - 1 );
}
y = u + ( Z >> b );

I know that this is not a strict mathematical proof, but it is confirmed by simulations.