# best implementation of a real-time, fixed-point iir filter with constant coefficients

I have implemented IIR filter in direct form, parallel form, and cascade form. The input in all the cases is a kronecker delta function. How to judge which implementation is better is in real time and why? What are the factors to be considered while selecting the best implementation? Please help.

• fixed point? floating point? word width? constant or time-varying coefficients? memory usage? computation time? how important to your design are these issues? Commented Feb 27, 2015 at 15:11
• i want in fixed point,the q formats are 8,16,24,the computation time should be very less,and the coeffficents are constant they are not time varying- Commented Feb 27, 2015 at 19:06
• dunno why someone else downvoted this. might have been light on fundamentals, but he provided enough information for the beginning of an answer. Commented Feb 28, 2015 at 23:45

what is fraction saving? can you write a code. so that i can understand more clearly?

Let's call the quantizer operator $$\operatorname{Quant}\{\cdot\}$$ . So the output of the quantizer, with $$v[n]$$ going in, is

$$y[n] = \operatorname{Quant}\{ v[n] \}$$

which we shall model as an additive error source:

$$y[n] = v[n] + q[n]$$

No matter how the quantizer works, we can always derive $$q[n]$$ in the run-time system with

$$q[n] \triangleq y[n] - v[n] = \operatorname{Quant}\{ v[n] \} - v[n]$$

Even though $$q[n]$$ is, strictly speaking, a deterministic function of $$v[n]$$, often, if the signal swings are much larger than the step size, $$\Delta$$, of the quantizer, we model $$q[n]$$ as some kinda random number and a "noise source". It's when the signal gets very small, on par with $$\Delta$$, that this assumption is no good, and that's when you need to consider adding dither before the quantization operation. I doubt you will need dither for your situation, but "ya never know".

The simplest quantizer to implement in a fixed-point 2's complement system is simply dropping the less significant bits that are to the right of the rounding point or "binary point":

$$\operatorname{Quant}\{ v[n] \} \triangleq \Delta \left\lfloor \frac{v[n]}{\Delta} \right\rfloor$$

$$\lfloor \cdot \rfloor$$ is the floor function, the same as the what the C function floor(): returning the most-positive integer that does not exceed the argument to the function. $$\lfloor u \rfloor$$ is the sole integer that satisfies:

$$u - 1 < \lfloor u \rfloor \le u$$

or

$$\lfloor u \rfloor \le u < \lfloor u \rfloor + 1$$

If $$u \ge 0$$, often we say that $$\lfloor u \rfloor$$ is the "integer part" of the real value $$u$$. The "fractional part" of $$u$$ is $$u - \lfloor u \rfloor$$. It is clear that the fractional part of $$u$$ is

$$0 \le u - \lfloor u \rfloor < 1$$

So the quantization error is the fractional part remaining after quantization is

$$q[n] = \Delta \left\lfloor \frac{v[n]}{\Delta} \right\rfloor - v[n]$$

That is the additive noise signal in terms of the input to the quantizer.

$$-1 < \frac{q[n]}{\Delta} = \left\lfloor \frac{v[n]}{\Delta} \right\rfloor - \frac{v[n]}{\Delta} \le 0$$

So we know that the additive quantization error, $$q[n]$$, has a range of

$$-\Delta < q[n] \le 0$$

Now, we assume (that's the only "rule-of-thumb") that essentially the value of $$q[n]$$ is equally likely in that range so then $$q[n]$$ is viewed as a random number with uniform p.d.f.:

$$p_q(\alpha) = \frac{1}{\Delta} \operatorname{rect}\left(\frac{\alpha}{\Delta} + \frac{1}{2}\right)$$

where $$\operatorname{rect}(\cdot)$$ is the standard unit rectangular function:

$$\operatorname{rect}(u) \triangleq \begin{cases} 1, & \text{if } \ -\frac{1}{2} < u \le \frac{1}{2} \\ 0, & \text{otherwise} \end{cases}$$

The expectation value or mean of $$q[n]$$ is $$-\frac{\Delta}{2}$$ (a half-step bias toward the negative direction) and the variance is $$\frac{\Delta^2}{12}$$. This is well known about uniform quantization.

In the frequency domain:

$$Y(z) = V(z) + Q(z)$$

(for frequency response, we substitute $$z = e^{j 2 \pi f T}$$ where $$T$$ is the sampling period and in reciprocal units as $$f$$. Normalized frequency (using the DTFT or "discrete-time Fourier transform") is $$\omega \triangleq 2 \pi f T$$ and ranges $$-\pi < \omega < \pi$$. In the DTFT, if the quantization error really is random, $$Q(e^{j \omega})$$ would have a Dirac impulse at $$\omega=0$$ of magnitude $$-\frac{\Delta}{2}$$ (to represent the DC bias of $$-\frac{1}{2}$$ LSB because we're always rounding down) and otherwise a flat spectrum from $$-\pi$$ to $$\pi$$, that would, because of Parseval's theorem, have area of $$\frac{\Delta^2}{12}$$. Maybe there's a $$2\pi$$ factor in that (which is why i don't usually like angular frequency over "ordinary frequency").

We know that:

$$\int\limits_{-\infty}^{+\infty} p_q(\alpha) \ d\alpha \ = \ 1$$

The mean of quantization error (DC component):

\begin{align} \mu_q & \triangleq \overline{q[n]} \\ & \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} q[n] \\ & = \int\limits_{-\infty}^{+\infty} p_q(\alpha) \ \alpha \ d\alpha \\ & = Q(e^{j 0}) \\ & = -\frac{\Delta}{2} \end{align}

The mean square of quantization error:

\begin{align} \overline{\left| q[n] \right|^2} & \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} |q[n]|^2 \\ & = \int\limits_{-\infty}^{+\infty} p_q(\alpha) \ \alpha^2 \ d\alpha \\ & = \frac{1}{2\pi}\int\limits_{-\pi}^{+\pi} \left| Q(e^{j \omega}) \right|^2 d\omega \\ & = \frac{\Delta^2}{3} \end{align}

The variance of quantization error (AC power):

\begin{align} \sigma_q^2 & \triangleq \overline{\left| q[n] - \overline{q[n]} \right|^2} \\ & \triangleq \lim_{N \to +\infty}\frac{1}{2N+1} \sum\limits_{n=-N}^{N} |q[n] - \overline{q[n]}|^2 \\ & = \int\limits_{-\infty}^{+\infty} p_q(\alpha) \ (\alpha-\mu_q)^2 \ d\alpha \\ & = \overline{\left| q[n] \right|^2} - \left|\overline{q[n]}\right|^2 \\ & = \frac{\Delta^2}{12} \end{align}

Now (to get back to "fraction saving") suppose $$v[n]$$ is derived from your input, $$x[n]$$, as so:

$$v[n] = x[n] - q[n-1]$$

So whatever error we added to $$v[n]$$ by quantizing it, we will subtract, before quantization, from the quantized sample immediately following.

(n.b.: Remember that $$q[n] \le 0$$ and when we "add" $$q[n]$$ in the quantization process, we are effectively subtracting whatever bits are to the right of the binary point by zeroing those bits, which is always rounding down. So when we "subtract" $$q[n-1]$$ from the quantizer input, $$x[n]$$, before quantization, we are adding those bits we dropped, extended on the left with zeroed bits, not sign extended.)

\begin{align} y[n] & = \operatorname{Quant}\{ v[n] \} \\ & = \operatorname{Quant}\{ x[n] - q[n-1] \} \\ & = \quad x[n] - q[n-1] + q[n] \\ & = \quad x[n] + (q[n] - q[n-1]) \\ \end{align}

or, in the frequency domain:

\begin{align} Y(z) & = X(z) \quad + \quad Q(z) - Q(z) z^{-1} \\ & = X(z) \quad + \quad (1 - z^{-1}) Q(z) \\ & = H_x(z) X(z) \quad + \quad H_q(z) Q(z) \\ \end{align}

So the transfer function from input $$x[n]$$ to output $$y[n]$$ is still just a wire:

$$H_x(z) = 1$$

But the transfer function from the additive quantization noise $$q[n]$$ to the output $$y[n]$$ is

$$H_q(z) = 1 - z^{-1} = \frac{z-1}{z}$$

So, like FIR filters, there is a pole at the origin $$z=0$$ (big hairy deel), but there is a zero right at the DC point, $$z=1$$. This kills any contribution of quantization error to the output at the frequency of 0. Even that $$-\frac{\Delta}{2}$$ bias, from always rounding down, is killed.

The magnitude of the noise to output frequency response is

\begin{align} |H_q(e^{j \omega})| & = \frac{|e^{j \omega} - 1|}{|e^{j \omega}|} \\ & = |e^{j \omega/2}| \frac{|e^{j \omega/2} - e^{-j \omega/2}|}{|e^{j \omega}|} \\ & = 2 \left| \sin\left(\frac{\omega}{2} \right) \right| \\ \end{align}

The transfer function, $$H_q(\cdot)$$, is 0 (or $$-\infty$$ dB) at DC ($$\omega=0$$) and when you're up at Nyquist (or $$\omega=\pi$$) the gain is 2 (or 6 dB). overall, there is no decrease in quantization noise energy, in fact, because energy is magnitude-squared, there is an increase in total quantization energy, but at frequencies below $$\frac{\pi}{3}$$, the quantization noise is reduced. For us audio guys, we like that because none of us hears super well above 8 kHz anyway.

But one thing it helps is that it destroys that limit cycle of IIR filters that gets them stuck on a non-zero DC value, even when the input truly goes to a dead zero.

Code for a biquad (Direct Form I, not Direct Form II which you don't wanna use, especially for fixed-point) that might look like that might be as below. To be clear in sign convention, the transfer function from input to output is:

$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 - a_1 z^{-1} - a_2 z^{-2}}$$

This corresponds to a DF1 difference equation as:

$$y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2]$$

With quantization it's:

$$y[n] = \operatorname{Quant}\{ b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] \}$$

With the noise-shaping (or fraction-saving) it's:

\begin{align} y[n] &= \operatorname{Quant}\{ b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] - q[n-1] \} \\ \\ &= \big( b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] - q[n-1] \big) + q[n] \\ \end{align}

or

$$q[n] = y[n] - \big(b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + a_1 y[n-1] + a_2 y[n-2] - q[n-1] \big)$$

Because, in terms of actual mathematical values, any stable IIR will have $$-2 < a_1 < +2$$ and $$-1 < a_2 < +1$$ ($$b_0$$, $$b_1$$, and $$b_2$$ are relative to each other so we worry less about their fixed range), so let's say that all our coefficients are all scaled up by a factor of $$2^{14} = 16384$$ (so that the equivalent integer value of the fixed-point number ranges from $$-2\cdot2^{14} = -32768 \ <$$ a1 $$= 2^{14} a_1 < \ +2\cdot2^{14} = +32768$$). So the quantization or rounding point is at bit 14 of the 32 bit number. The coefficients are also stored as 16-bit int16_t. We cast the coefficient to int32_t because we want coefficient times state to be an int32_t, even if both multiplier and multiplicand going in are both int16_t.

Version 1:

#include <stdint.h>

typedef struct
{
int16_t num_samples;            // number of samples in sample block (same for input and output)
int16_t *input, *output;        // input and output sample block pointers
int16_t x1, x2, y1, y2;         // Direct Form 1 2nd-order IIR states
int16_t saved_fraction;         // noise-shaping state
int16_t a1, a2, b0, b1, b2;     // 2nd-order IIR coefficients

{
int16_t x1 = this_biquad->x1;           // get states

int16_t b0 = this_biquad->b0;           // get coefficients

{
int16_t x = *input++;

accumulator += (int32_t)b0*x;           // an optimized compiler will figure this out
accumulator += (int32_t)b1*x1;          //  and perform only a 16x16=32bit multiply-accumulate
accumulator += (int32_t)b2*x2;
accumulator += (int32_t)a1*y1;
accumulator += (int32_t)a2*y2;

if (accumulator > 0x1FFFFFFFL)
{
accumulator = 0x1FFFFFFFL;          // clip value
}
if (accumulator < -0x20000000L)
{
accumulator = -0x20000000L;         // clip value
}

int16_t y =  (int16_t)(accumulator>>14);        // always rounding down

x2 = x1;                                // bump the states over
x1 = x;
y2 = y1;
y1 = y;

accumulator &= 0x00003FFFL;
// keep the fractional bits that you dropped for the next sample otherwise clear the accumulator

*output++ = y;
}

this_biquad->x1 = x1;                   // save states
}


Now another trick, to speed things up, is to adjust the a1 coefficient by subtracting the scaled fixed-point version of 1 (which, in this case, is $$2^{14}$$). Then the code might look like:

Version 2:

#include <stdint.h>

typedef struct
{
int16_t num_samples;            // number of samples in sample block (same for input and output)
int16_t *input, *output;        // input and output sample block pointers
int16_t x1, x2, y1, y2;         // Direct Form 1 2nd-order IIR states
int16_t saved_fraction;         // noise-shaping state
int16_t a1, a2, b0, b1, b2;     // 2nd-order IIR coefficients

{
int16_t x1 = this_biquad->x1;           // get states
int32_t accumulator = ((int32_t)y1<<14) + (int32_t)(this_biquad->saved_fraction);

int16_t b0 = this_biquad->b0;           // get coefficients
int16_t a1 = this_biquad->a1 - 0x4000;  // subtract 1 from a1 coefficient

{
int16_t x = *input++;

accumulator += (int32_t)b0*x;           // an optimized compiler will figure this out
accumulator += (int32_t)b1*x1;          //  and perform only a 16x16=32bit multiply-accumulate
accumulator += (int32_t)b2*x2;
accumulator += (int32_t)a1*y1;          // this coefficient is actually a1 - 1
accumulator += (int32_t)a2*y2;

if (accumulator > 0x1FFFFFFFL)
{
accumulator = 0x1FFFFFFFL;          // clip value
}
if (accumulator < -0x20000000L)
{
accumulator = -0x20000000L;         // clip value
}

int16_t y =  (int16_t)(accumulator>>14);        // always rounding down

x2 = x1;                                // bump the states over
x1 = x;
y2 = y1;
y1 = y;

*output++ = y;
}

this_biquad->x1 = x1;                   // save states
}


Note the missing bit-masking instruction at the bottom of the loop. So you have the left bits having an extra value of $$y[n-1]$$ in there. That's why we have to fix a1 by subtracting 1 (again scaled by $$2^{14}$$) from it. That does affect the range of a1 by bumping it over. Now it's $$-2 < a_1 - 1 < +2$$ or $$-1 < a_1 < +3$$ but the range of all of the other coefficients are unchanged.

Below is code that does two things:

1. Doubles the word size for everything. Now we got 28 fractional bits and our fixed-point coefficients can range from -8.00000000 to +7.99999999. And 24 dB of headroom.
2. Coefficients are integers scaled up $$2^{28}$$ from their mathematical value.
3. Cascades as many biquad sections (still Direct Form 1) as you can handle with memory and speed resources. Each section has their own saved_fraction state.

This is extremely high-quality processing. Lot'sa precision. Better than single-precision floating point.

#include <stdint.h>

/*
int32_t filter_states[3*num_sections+2];        // [x1, x2, saved_fraction, y1, y2]        for first section
// then                    [x1, x2, saved_fraction, y1, y2]
// for each following section.  Output states (y1,y2) for
// a given section become input states (x1,x2) in the following
// section.  Filter states should be set to zero to begin, but

int32_t filter_coefficients[5*num_sections];    // b0, b1, b2, a1, a2 for each section.  Coefficients are
//  integers scaled up by 2^28 from their mathematical value
//  that can range from -8.00000000 to +7.99999999 .

int32_t output[num_samples];                    // This processes samples in blocks of size num_samples.
int32_t input[num_samples];
*/

int32_t* filter_states,
int32_t* filter_coefficients,
int num_samples,
int32_t* output,
int32_t* input )
{
int32_t* state_ptr = filter_states;
int32_t* coef_ptr = filter_coefficients;

int32_t x1, x2, y1, y2;
int64_t accumulator;

int32_t* in = input;                    // this input pointer used only for first section

for (int k=num_sections; k>0; k--)
{
int32_t* out = output;

x1 = *state_ptr++;                  // get states
x2 = *state_ptr++;
accumulator = (int64_t)(*state_ptr++);
y1 = *state_ptr++;
y2 = *state_ptr;
state_ptr -= 4;                     // backup to point to x1 state

int32_t b0 = *coef_ptr++;           // get coefficients
int32_t b1 = *coef_ptr++;
int32_t b2 = *coef_ptr++;
int32_t a1 = *coef_ptr++;
int32_t a2 = *coef_ptr++;

for (int n=num_samples; n>0; n--)
{
int32_t x = *in++;

accumulator += (int64_t)b0*x;           // an optimized compiler will figure this out
accumulator += (int64_t)b1*x1;          //  and perform only a 32x32 = 64 bit multiply-accumulate
accumulator += (int64_t)b2*x2;
accumulator += (int64_t)a1*y1;
accumulator += (int64_t)a2*y2;

if (accumulator > 0x07FFFFFFFFFFFFFFLL)
{
accumulator = 0x07FFFFFFFFFFFFFFLL;      // clip value
}
if (accumulator < -0x0800000000000000LL)
{
accumulator = -0x0800000000000000LL;     // clip value
}

int32_t y =  (int32_t)(accumulator>>28);     // point of quantization, always rounding down

x2 = x1;                                     // bump the states over
x1 = x;
y2 = y1;
y1 = y;

accumulator &= 0x000000000FFFFFFFLL;     // keep the fractional bits that were dropped for
// the next sample, otherwise clear the accumulator
*out++ = y;
}

*state_ptr++ = x1;               // save states
*state_ptr++ = x2;
*state_ptr++ = (int32_t)accumulator;

in = output;                     // subsequent sections use the same output sample block as input
}

*state_ptr++ = y1;                   // save last section output states
*state_ptr = y2;
}

• One of the most useful answers ever ;) Thanks Robert!
– jojeck
Commented Mar 12, 2015 at 9:32
• tried to hit the theory and the code to fully explain it, @jojek . most noise shaping is cheap, easy, and effective. fraction saving has a particularly high payoff vs cost ratio. it costs almost nothing and it totally kills any quantization error at DC, which totally kills that common limit-cycle where the output of your filter gets stuck on a non-zero DC value after the input to the same filter goes to a dead zero. Commented Mar 12, 2015 at 14:12
• contact me by email (my address is easy to find if you google me), and we'll work out a fee and arrange a consulting contract. it's pretty clear in C language and you'll have to figure out how to do integer or fixed-point (and perhaps logical) operations in MATLAB. actually, you can do this without fixed-point, and to implement the quantizer, you need floor( ) in MATLAB. then subtract to get the quantization error. Commented Mar 26, 2015 at 18:18
• well, the area of the curve under $\left|Q\left( e^{j\omega} \right)\right|^2$ for $0 < \omega < +\pi$ (i want to exclude the DC component due to always rounding down) is $\frac{\Delta^2}{24}$ or maybe i am missing a factor of $2 \pi$ in that because of the use of angular frequency. so the height of that is either $\frac{\Delta^2}{24\pi}$ or $\frac{\Delta^2}{12}$. (i have to look up exactly how power spectrum is done for angular frequency.) Commented Mar 31, 2015 at 19:21
• Hi robert, that's an excellent answer! I have a minor doubt that you said at frequencies below $\pi/2$, the quantization noise is reduced. According to the frequency response $2|\sin(\omega/2)|$ I think it should be $\pi/3$, am I right? Commented Aug 16, 2021 at 4:37

i want in fixed point,the q formats are 8,16,24,the computation time should be very less,and the coeffficents are constant they are not time varying.

i guess i would recommend, Direct Form I. use an accumulator that with word width equal to the sum of word widths of signal and coefficient. so there's no accumulation of quantization error until the last step of computing the output signal word. in that last step where you round or quantize from the accumulator to the output signal, i would suggest what is commonly called "fraction saving", but could be called "first-order noise shaping with round-down quantizer": whatever bits you lose (by truncation) in rounding down, add those zero-extended bits into the accumulator in the following sample. this way your rounding error at DC is zero. infinite S/N at DC, lowered S/N for low frequencies and slightly higher S/N for frequencies above Nyquist/3 . and it kills a certain limit cycle that can be an annoyance for recursive fixed-point filters.

and it's as cheap, computationally, as any other form (and cheaper than most other forms). it does require 2N+2 states for an order 2N string of cascaded biquad sections. that's 2 more than the canonical Direct Form II.

• what is fraction saving.can you write a code.so that i can understand more clearly Commented Feb 27, 2015 at 20:12