Single pole IIR filter, fixed point design

Question

We want to do a fixed point implementation of the single pole IIR filter:

$y[n] = a\cdot x[n] + (1 - a)\cdot y[n-1] \quad ;\qquad 0<a<1$

What are the main design considerations to think about ? In particular, is the fixed point design more challenging for smaller a or larger a ?

Answer 1

One thing to consider when implementing an IIR filter, whatever the order, is quantization and limit cycles.

Let me show you with a quick example with your original filter

$y[n] = a*x[n]+(1-a)*y[n-1]$

Let a = 0.005 and say that we use 16-bit signed coefficients.

$a_{fixedpoint} = a * 32768 = 164$

Let's assume that the input and output are 16-bit integers. They will go from -32768 to 32767. The previous output $y[n-1]$ is also a 16-bit integer in this scenario. The implementation would look like this, with temp being a 32-bit integer. All other variables are 16-bit integers.

$ temp = a_{fixedpoint}*x[n] + (32768 - a_{fixedpoint})*y[n-1]$

$y[n] = round(temp/32768)$

Let's try this implementation with a step response

It looks good, right? Perhaps we should zoom in.

We don't reach 32767 even though the gain of the original filter is 1. This is called a limit cycle. For order-1 IIR filters, limit cycles take the form of a DC offset that you can't get rid of. For IIR filters with complex poles, limit cycles take the form of a small oscillation that you can't get rid of. The problem is caused by the fact that $y[n-1]$ is stored on 16 bits. If we increase the resolution of $y[n-1]$ to 32 bits, it would solve the problem. Another solution would be to use fraction-saving.

If the input would go from 32767 to 0, we would have the same problem, the output would get close to 0 but would not actually reach 0.

Increased resolution

If instead of storing the output on 16 bits, we store it on 24 bits and name this variable acc (increased resolution of y[n]).

$ temp = a_{fixedpoint}*x[n]*256 + (32768 - a_{fixedpoint})*acc$

$ acc = temp/32768$

$y[n] = round(acc/(256))$

As we can see, with the increased resolution we can reach 32767 even though the coefficients are still on 16 bits. However the price to pay is having intermediate variables with greater resolution, the calculations might be slower.

fraction-saving

Finally there's a clever approach to minimize the added resolution needed for an IIR, especially when the poles are close to the unit circle. It is called fraction-saving. The idea is to memorize the round-off error of the previous sample and apply it to the next calculation to reduce the quantization effects. An additional variable is needed called error, but you don't need to increase the resolution of the previous output. So if x is represented on 16 bits, then y[n] can be stored on 16 bits too. The error variable would be represented with 16 bits too in this example. As Robert Bristow-Johnson pointed out, fraction-saving is a form of noise-shaping

$ temp = a_{fixedpoint}*x[n] + (32768 - a_{fixedpoint})*y[n-1] - error[n-1]$

$y[n] = round(temp/32768)$

$ error[n] = y[n]*32768 - temp$

Answer 2

this is my original posting about implementing a fixed-point DC blocking first-order HPF using fraction saving.

Lyons and Yates later did an IEEE Sig Proc article about DC blocking filters where this was one of the topics.

// let's say sizeof(short) = 2 (16 bits) and sizeof(long) = 4 (32 bits)
short x[], y[];
long acc, A, prev_x, prev_y;
double pole;
unsigned long n, num_samples;
pole = 0.9999;
A = (long)(32768.0*(1.0 - pole));
acc = 0;
prev_x = 0;
prev_y = 0;
for (n=0; n<num_samples; n++)
    {
    acc   -= prev_x;
    prev_x = (long)x[n]<<15;
    acc   += prev_x;
    acc   -= A*prev_y;
    prev_y = acc>>15;               // quantization happens here
    y[n]   = (short)prev_y;
    // acc has y[n] in upper 17 bits and -e[n] in lower 15 bits
    }

Answer 3

In particular, is the fixed point design more challenging for smaller a or larger a ?

Smaller is worse: the closer $a$ gets to $0$, the closer the pol movess to the unit circle and the more time domain ringing you have.

What are the main design considerations to think about ?

A first order low-pass IIR filter is relatively benign. The sum absolute sum of the impulse response is unity so you can't clip the output. One thing to watch out for are potentially limit cycles. This can be controlled by the way you round: "round towards zero" is the safest option.

Answer 4

It used to be that multiplication was a lot more expensive than addition or bitwise operations. For integer (and fixed point) implementation, this can be taken exploited like this:

When $a=1/2$ an extremely efficient implementation can be made with a single add and a single shift.

$$ y[n] = (x[n] + y[n-1]) >> 1 $$

Similarly, fractions with a power of two in the denominator can also be done efficiently. For instance, $a = 7/8$ can be done as:

$$ y[n] = ((x[n]<<3) - x[n] + y[n-1]) >> 3 $$

There are many other possibilities as well. I tend to go with the $a=1/2$ option a lot.

Obviously, how much headroom you have is a major consideration, so doing something like $a = 255/256$ would need 8 bits of headroom, whereas having fixed point fractions doesn't suffer from that, but it costs you a multiplication or two. Note the equation can also be rewritten as:

$$ y[n] = y[n-1] + a * ( x[n] - y[n-1] ) $$

This saves a multiplication at the cost of a subtraction.

Answer 5

Big simplification can be done on equations of first answer (fraction-saving) as follows:

Original:

temp = afixedpoint∗x[n] + (32768-afixedpoint)∗y[n-1] - error[n-1]
y[n] = round(temp/32768)
error[n] = y[n]∗32768−temp

Proposed:

temp = temp + afixedpoint*(x[n]-y[n-1])
y[n] = round(temp/32768)

Coded was tested giving same results

Proof

The (32768-a) term in original equation can be simplified with this:

a*x+(k-a)*y = a*x+k*y-a*y = a*(x-y)+k*y (with k = 32768)

Using it in the original equation:

temp = a*(x[n]-y[n-1]) + k*y[n-1] - e[n-1]
y[n] = round(temp/k)
e[n] = k*y[n]-temp

Substituting error term:

temp = a*(x[n]-y[n-1]) + k*y[n-1] - (k*y[n-1]-temp[n-1])
temp = a*(x[n]-y[n-1]) + k*y[n-1] - k*y[n-1] + temp[n-1]
temp = a*(x[n]-y[n-1]) + temp[n-1]

Simplifying:

temp = temp + a*(x[n]-y[n-1])
y[n] = round(temp/k)