Aversive/Modules/control system/filters/biquad

Description

This filter is a general IIR digital filter. All IIR and FIR filters can be implemented using this module. The general equation of this filter (in Z transform) is :

        b0/do + b1/do*z^-1 + b2/do*z^-2                   1                       1              
H(z) = ---------------------------------      with do = -----------------  dr = -----------------
         1    + a1/dr*z^-1 + a2/dr*z^-2                  2^out_shift             2^reac_shift     

The overall gain is not adjustable, for simplicity reasons. This implementation, in direct form I is probably a good compromise between memory, execution time and immunity to quantification errors, and so on. The choice to stay in fixed point has been done for a lot of reasons, the first one is better precision.

For more details, you can consult this page :[1]

use as derivation, integration, or just as a simple coefficient

The biquad can be used as general derivation or integration block.

Derivation example :

struct biquad_filter derivation;
// derivation 
biquad_init(&derivation);
biquad_set_numerator_coeffs(&derivation, 1,-1,0); // this is a discrete derivation : O(n) = I(n) - I(n-1)
// no need to initialize denominator coeffs to 0, init has done it

Integration is done this way :

struct biquad_filter integration;
// integration
biquad_init(&integration);
biquad_set_denominator_coeffs(&integration, 1,0); // this is a discrete integration : O(n) = I(n) - O(n-1)
// no need to initialize numerator coeffs to 0, init has done it

Coefficient :

struct biquad_filter coeff;
// multiplication
biquad_init(&coeff);
biquad_set_numerator_coeffs(&derivation, 5,0,0); // multiply by 5
biquad_set_divisor_shifts(&derivation, 0, 6);    // and divide by 2^6 = 64

example of a filter design procedure

TODO : add example

First you need to calculate coefficients, and choose filter order.

For this you can use any filter design method.

• pole qnd zero placement : [2]
• simple equations for second order : Biquad : [3]
• chebyshev, butterworth, ... filter design methods for higher orders: you will need a design tool. you can use [4] or in matlab [5] You will also have to splitt your filter into second order biquads, if you have a higher order.

Once you have the filter coefficients, you can determine the value of divisor_shift. This shift should be as big as possible to have enough precision on the coefficients. This is determined by your biggest coefficient.

Concrete example :

Sampling rate is 1000 Hz. We want a LPF(low pass) with 50 Hz cutoff. We will use a second order butterworth

With these parameters, [6] gives back :

Recurrence relation:
y[n] = (  1 * x[n- 2])
     + (  2 * x[n- 1])
     + (  1 * x[n- 0])

     + ( -0.6413515381 * y[n- 2])
     + (  1.5610180758 * y[n- 1])

which means (beware on the sign on y terms!):

b2' = 1
b1' = 2
b0' = 1
a2' = 0.6413515381
a1' = -1.5610180758

Our max coeff is 2, and the max span of the s16 coeffs is +-32767. We have now to choose do and dr. For the b coefficients, it's relatively easy, because we have round coeffs, so we will chose do as low as possible to avoid overflow : do = 2^0 For the a coefficients, we do not have round values, so we will have to choose a relatively high dr, in order to have precision. But this is a tradeoff with the overflow capabilities. We choose a mid value of dr = 2^8. any shift between 3 and 15 can be used.

this gives :

dr = 2^8, do = 2^0
b2 = 1*2^0        =   1      >>  rounded b2' = 1
b1 = 2*2^0        =   2      >>  rounded b1' = 2
b0 = 1*2^0        =   1      >>  rounded b0' = 1
a2 =  0.6...*2^8  =~  164    >>  rounded a2' = 164/2^8 = 0.640625
a1 =-1.56...*2^8  =~  -400   >>  rounded a1' =-400/2^8 =-1.5625

With this method we can minimize the quantification error on the coeffs by increasing dr.

One problem is remaining : this filter has a gain of 50 : gain at dc  : mag = 4.979245121e+01

For getting rid of this, we will change b0,b1, and b2 (and do):

b2' = 1/50
b1' = 2/50
b0' = 1/50

for implementing this, we will have to make do higher. We choose do= 2^10

dr = 2^8, do = 2^8
b2 = 0.02*2^8   =~ 5      >>  rounded b2' = 5/2^8 = 0.019531
b1 = 0.04*2^8   =~ 10     >>  rounded b1' = 10/2^8= 0.039063
b0 = 0.02*2^8   =~ 5      >>  rounded b0' = 5/2^8 = 0.019531
a2 =  0.6...*2^8  =~  164    >>  rounded a2' = 164/2^8 = 0.640625
a1 =-1.56...*2^8  =~  -400   >>  rounded a1' =-400/2^8 =-1.5625

Finally we write the following code :

biquad_init(&filter);
biquad_set_numerator_coeffs  (&filter, 5, 10, 5);
biquad_set_deniminator_coeffs(&filter,-400, 164);
biquad_set_divisor_shifts(&filter, 8, 8);
...
output = biquad_do_filter(& filter1, input);

The results of this filter can be shown graphically:

in the spectrum, this gives following :

Please beware of overflows !! There is no security in this filter (actually), and even if the input is a s32, you should use the minimal necessary input (generally the value of your input samples from ADC, 10 bits) and not multiplicate this before the filter!

You can also make tradeoff in reducing do and/or dr, but you will loose precision on the coeff rounding.

Have a happy filter design !