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

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
} biquad;

void biquad_process_samples(biquad* this_biquad)
{
    int16_t x1 = this_biquad->x1;           // get states
    int16_t x2 = this_biquad->x2;
    int16_t y1 = this_biquad->y1;
    int16_t y2 = this_biquad->y2;
    int32_t accumulator = (int32_t)(this_biquad->saved_fraction);

    int16_t b0 = this_biquad->b0;           // get coefficients
    int16_t b1 = this_biquad->b1;
    int16_t b2 = this_biquad->b2;
    int16_t a1 = this_biquad->a1;
    int16_t a2 = this_biquad->a2;

    int16_t* input = this_biquad->input;
    int16_t* output = this_biquad->output;

    for (int16_t n=this_biquad->num_samples; n>0; n--)
        {
        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
    this_biquad->x2 = x2;
    this_biquad->y1 = y1;
    this_biquad->y2 = y2;
    this_biquad->saved_fraction = (int16_t)accumulator;
}

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
} biquad;

void biquad_process_samples(biquad* this_biquad)
{
    int16_t x1 = this_biquad->x1;           // get states
    int16_t x2 = this_biquad->x2;
    int16_t y1 = this_biquad->y1;
    int16_t y2 = this_biquad->y2;
    int32_t accumulator = ((int32_t)y1<<14) + (int32_t)(this_biquad->saved_fraction);

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

    int16_t* input = this_biquad->input;
    int16_t* output = this_biquad->output;

    for (int16_t n=this_biquad->num_samples; n>0; n--)
        {
        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
    this_biquad->x2 = x2;
    this_biquad->y1 = y1;
    this_biquad->y2 = y2;
    this_biquad->saved_fraction = (int16_t)accumulator & 0x3FFF;
}

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 
                                                // must be left unchanged between calls to cascaded_biquads().

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];
*/

void cascaded_biquads( int num_sections, 
                       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;
}

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.