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

Question

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.

Answer 1

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] \} $$

I am sorta assuming that short are 16-bit words and long are 32-bit words. 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 short, but we are casting them to long upon loading so that the cast is not continuously necessary in the inner loop. We want coefficient times state to be a long, even if both multiplier and multiplicand going in are short.

Version 1:

long accumulator;

long a1, a2, b0, b1, b2;   /* 2nd-order IIR coefficients */

short state_x1, state_x2, state_y1, state_y2, state_error;

short x, y;

short get_input(void);

void put_output(short y);

short i, num_samples;
         /* num_samples is the number of samples per block */

/*  load up your coefficients, cast all to long */

/*  load up your states from whereever */

accumulator = (long)state_error;

for (i=0; i<num_samples; i++)
    {
    x = get_input();

    accumulator += b0*x;
    accumulator += b1*state_x1;
    accumulator += b2*state_x2;
    accumulator += a1*state_y1;
    accumulator += a2*state_y2;

    if (accumulator > 0x1FFFFFFF)
       {
       accumulator = 0x1FFFFFFF;     /* clip value */
       }
    if (accumulator < -0x20000000)
       {
       accumulator = -0x20000000;    /* clip value */
       }

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

    state_x2 = state_x1;             /* bump the states over */
    state_x1 = x;
    state_y2 = state_y1;
    state_y1 = y;

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

    put_output(y);
    }

state_error = (short)accumulator;

/*  save your states back to whereever  */

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:

long accumulator;

long a1, a2, b0, b1, b2;   /* 2nd-order IIR coefficients */

short state_x1, state_x2, state_y1, state_y2, state_error;

short x, y;

short get_input(void);

void put_output(short y);

short i, num_samples;
         /* num_samples is the number of samples per block */

/*  load up your coefficients, cast all to long */

/*  load up your states from whereever */

a1 = a1 - 0x00004000;    /* adjust a1 coefficient */
/* 0x00004000 is "1" in our scaled-up fixed-point coefficients */

accumulator = (long)state_error + (((long)state_y1)<<14);

for (i=0; i<num_samples; i++)
    {
    x = get_input();

    accumulator += b0*x;
    accumulator += b1*state_x1;
    accumulator += b2*state_x2;
    accumulator += a1*state_y1;
    accumulator += a2*state_y2;

    if (accumulator > 0x1FFFFFFF)
       {
       accumulator = 0x1FFFFFFF;     /* clip value */
       }
    if (accumulator < -0x20000000)
       {
       accumulator = -0x20000000;    /* clip value */
       }

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

    state_x2 = state_x1;             /* bump the states over */
    state_x1 = x;
    state_y2 = state_y1;
    state_y1 = y;

    put_output(y);
    }

state_error = (short)(accumulator & 0x00003FFF);

/*  save your states back to whereever  */

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.

Answer 2

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.