0
$\begingroup$

I've been trying to make a digital filter. I've used the following transfer function to get my digital coefficients https://gyazo.com/00195ef119f0b2ff93fa5855715abdc9

Which i've then transfered to code using the direct form 1 with 2 2nd order sections. Which you can see down below in the code.

To get this into smaller coefficients i first divided each section by a10 and a20 to get smaller coefficients and making a10 and a20 1. The undivided coefficients can still be seen in the code as comments.

To then get the code in fixed point i tried multiplying each coefficient by 2^15 to use 16 bit integers and then dividing my end result by the same amount. However this does not seem to work.

Voltage is the reading of the ADC which is the input of the filter. And this is a 4th order lowpass butterworth filter with a -3dB at 500hz

And the delays in the for loop can probably be done in a different faster way but i'm not sure how else to do this.

Also not sure if this is important but this is done on a 8bit atmega chip.

TLDR: How do i implement this in fixed point and can the for loop for the delays be done faster?

Kind regards Chris

double xA[3] = {0};
double yA[3] = {0};
double xB[3] = {0};
double yB[3] = {0};
double y2 = 0;
double y = 0;
xA[0] = voltage;
xB[0] = voltage;
    //2nd order

/*  
double a10 = 1.842071000613218e+09  
double a11 = -3.180260791197821e+09             
double a12 = 1.377668208188961e+09
double b10 = 9.869604401089357e+06
double b11 = 1.973920880217871e+07
double b12 = 9.869604401089357e+06

double a20 = 1.706052605083395e+09
double a21 = -3.180260791197821e+09
double a22 = 1.513686603718784e+09
double b20 = 9.869604401089357e+06
double b21 = 1.973920880217871e+07
double b22 = 9.869604401089357e+06


 */

// Coefficients for 2 2nd order filters.       
double a10 = 1      
double a11 = -1.72645939822                 
double a12 = 0.74789093782
double b10 = 0.0053578849
double b11 = 0.0107157698
double b12 = 0.0053578849

double a20 = 1
double a21 = -1.86410476542
double a22 = 0.88724497662
double b20 = 0.0057850528   
double b21 = 0.0115701056
double b22 = 0.0057850528



    y2 = b10* xA[0]+ b11 *xA[1]+ b12*xA[2]- a11*yA[1] - a12*yA[2]; // 1st 2nd order
    y2 = y2/a10;
    yA[0] = y2;

    y = b20* xB[0]+ b21 *xB[1]+ b22*xB[2]- a21*yB[1] - a22*yB[2];                       //2nd 2nd order
    y = y/a20;
    yB[0] = y;


    for(uint8_t i = 2; i > 0; i--)                             // delay for 1st 2nd order
        xA[i] = xA[i-1];

        for(uint8_t i = 2; i > 0; i--)
        yA[i] = yA[i-1];

        for(uint8_t i = 2; i > 0; i--)                        // delay for 2nd 2nd order
        xB[i] = xB[i-1];

        for(uint8_t i = 2; i > 0; i--)
        yB[i] = yB[i-1];

      voltage = y;

My fixed point attempt.

double xA[3] = {0};
double yA[3] = {0};
double xB[3] = {0};
double yB[3] = {0};
int32_t y2 = 0;
int32_t y = 0;
xA[0] = voltage;
xB[0] = voltage;

// Coefficients for 2 2nd order filters. 
// Times 2^15 for fixed point   
int16_t a10 = 16384

int16_t a11 = -28286                    
int16_t a12 = 12253
int16_t b10 = 88
int16_t b11 = 176
int16_t b12 = 88

int16_t a20 = 16384
int16_t a21 = -30541
int16_t a22 = 14537
int16_t b20 = 95
int16_t b21 = 190
int16_t b22 = 95



    y2 = b10* xA[0]+ b11 *xA[1]+ b12*xA[2]- a11*yA[1] - a12*yA[2]; // 1st 2nd order
    y2 = y2/a10;
    yA[0] = y2;

    y = b20* xB[0]+ b21 *xB[1]+ b22*xB[2]- a21*yB[1] - a22*yB[2];                       //2nd 2nd order
    y = y/a20;
    yB[0] = y;


    for(uint8_t i = 2; i > 0; i--)                             // delay for 1st 2nd order
        xA[i] = xA[i-1];

        for(uint8_t i = 2; i > 0; i--)
        yA[i] = yA[i-1];

        for(uint8_t i = 2; i > 0; i--)                        // delay for 2nd 2nd order
        xB[i] = xB[i-1];

        for(uint8_t i = 2; i > 0; i--)
        yB[i] = yB[i-1];

      voltage = y*0.00006103515;// divided by 16384 or 2^14

3rd edit with ur example. unsure about these lines in the example that you gave. and i removed the forloop at the bottom. The entire process is done using an interupt timer thus not using the forloop.

state_error = (short)(accumulator & 0x00003FFF); accumulator = (long)state_error + (((long)state_y1)<<14);

Also the adc value cast to voltage should probably done at the bottom and instead the res from adc read should be used for the calculations as it makes no sense to use the voltage for calculations to me.

//ADC_voltage:
#define MAX_VALUE_ADC   2047                               // only 11 bits are used
#define VCC         3.30
#define VREF        (((double) VCC) / 1.6)                 // maximal value is 2.06125 V
//DAC
#define MAX_VALUE_DAC 4095
#define CAL_DAC 1.000  // Calibration value DAC

double xA[3] = {0};
double yA[3] = {0};
double xB[3] = {0};
double yB[3] = {0};
int32_t y2 = 0;
int32_t y = 0;

// Times 2^14 for fixed point   
int32_t a10 = 16384

int32_t a11 = -28286                    
int32_t a12 = 12253
int32_t b10 = 88
int32_t b11 = 176
int32_t b12 = 88

int32_t a20 = 16384
int32_t a21 = -30541
int32_t a22 = 14537
int32_t b20 = 95
int32_t b21 = 190
int32_t b22 = 95


ISR(TCE0_OVF_vect)              //TIMER Interupt
{

    double  voltage ;       // contains read in voltage (mVolts)
    int16_t res;            // contains read in adc value
    int16_t BinaryValue;    // contains value to write to DAC


    res = read_adc();           //Read out the ADC (PIN A2)
    voltage = (double) res * 1000 * VREF / (MAX_VALUE_ADC + 1);                 // Measured voltage in Volts.

    y2 = b10* xA[0]+ b11 *xA[1]+ b12*xA[2]- a11*yA[1] - a12*yA[2]; // 1st 2nd order
    y2 = y2/a10;    // is this necessary? without fixed point its 1/1
    yA[0] = y2;
    if (y2 > 0x1FFFFFFF)
       {
       y2 = 0x1FFFFFFF;     /* clip value */
       }
    if (y2 < -0x20000000)
       {
       y2 = -0x20000000;    /* clip value */
       }

    y = b20* xB[0]+ b21 *xB[1]+ b22*xB[2]- a21*yB[1] - a22*yB[2];                       //2nd 2nd order
    y = y/a20;    // same for this one
    yB[0] = y;
    if (y > 0x1FFFFFFF)
       {
       y = 0x1FFFFFFF;     /* clip value */
       }
    if (y < -0x20000000)
       {
       y = -0x20000000;    /* clip value */
       }

    xA[0] = voltage;                                                     // 1st 2nd order delays
    xA[2] = xA[1];
    xA[1] = xA[0];

    yA[2] = yA[1];
    yA[1] = yA[0];


    xB[0] = voltage;                                                    // 2d 2nd order delays
    xB[2] = xB[1];
    xB[1] = xB[0];

    yB[2] = yB[1];
    yB[1] = yB[0];

    voltage = (short)(y >>14);                                          // divided by 16384 or 2^14     voltage = y>>14
    BinaryValue =  voltage*((MAX_VALUE_DAC)/(VCC))*0.001*CAL_DAC ;      // Bitvalue
    DACB.CH0DATA = BinaryValue ;                                    //write &USBDataIn to DAC (PIN A10)
    while (!DACB.STATUS & DAC_CH0DRE_bm);

I've changed the states to ints and now do the actual math with the adc res value instead of the voltage value which probably makes more sense. I've also removed the clipping if's. At this point i don't know how to improve any more.

//ADC_voltage:
#define MAX_VALUE_ADC   2047                               // only 11 bits are used
#define VCC         3.30
#define VREF        (((double) VCC) / 1.6)                 // maximal value is 2.06125 V
//DAC
#define MAX_VALUE_DAC 4095
#define CAL_DAC 1.000  // Calibration value DAC

int32_t xA[3] = {0};
int32_t yA[3] = {0};
int32_t xB[3] = {0};
int32_t yB[3] = {0};
int32_t y2 = 0;
int32_t y = 0;

// Times 2^14 for fixed point   
int32_t a10 = 16384
int32_t a11 = -28286                    
int32_t a12 = 12253
int32_t b10 = 88
int32_t b11 = 176
int32_t b12 = 88

int32_t a20 = 16384
int32_t a21 = -30541
int32_t a22 = 14537
int32_t b20 = 95
int32_t b21 = 190
int32_t b22 = 95

ISR(TCE0_OVF_vect)              //TIMER Interupt
{
    double  voltage ;       // contains read in voltage (mVolts)
    int16_t res;            // contains read in adc value
    int16_t BinaryValue;    // contains value to write to DAC
    res = read_adc();           //Read out the ADC (PIN A2)
    voltage = (double) res * 1000 * VREF / (MAX_VALUE_ADC + 1);     // Measured voltage in Volts.

    y2 = b10* xA[0]+ b11 *xA[1]+ b12*xA[2]- a11*yA[1] - a12*yA[2]; // 1st 2nd order
    y2 = y2/a10;    // is this necessary? without fixed point its 1/1
    yA[0] = y2;

    y = b20* xB[0]+ b21 *xB[1]+ b22*xB[2]- a21*yB[1] - a22*yB[2];                       //2nd 2nd order
    y = y/a20;    // same for this one
    yB[0] = y;

    xA[0] = res;                                                     // 1st 2nd order delays
    xA[2] = xA[1];
    xA[1] = xA[0];

    yA[2] = yA[1];
    yA[1] = yA[0];

    xB[0] = res;                                                    // 2d 2nd order delays
    xB[2] = xB[1];
    xB[1] = xB[0];

    yB[2] = yB[1];
    yB[1] = yB[0];

    res = (short)(y >>14);                                          // divided by 16384 or 2^14     voltage = y>>14
    voltage = (double) res * 1000 * VREF / (MAX_VALUE_ADC + 1);     // Measured voltage in Volts.

    BinaryValue =  voltage*((MAX_VALUE_DAC)/(VCC))*0.001*CAL_DAC ;      // Bitvalue
    DACB.CH0DATA = BinaryValue ;                                    //write &USBDataIn to DAC (PIN A10)
    while (!DACB.STATUS & DAC_CH0DRE_bm);
$\endgroup$
  • $\begingroup$ Why does your integer code not seem to work? Be specific. $\endgroup$ – hotpaw2 May 17 at 15:14
  • $\begingroup$ I'm assuming overflow issues from what i've read online. and i'm unsure if multiplying the a10 to b22 by 2^15 is the right way to do it. as this is just one of the ways to do it that i've found mentioned online. Regards Chris - Edit The coefficients a11 and a21 don't even fit in a 16 bit int if i multiply by 2^15. Which i'm assuming is part of the problem $\endgroup$ – Chriserke May 17 at 15:31
  • $\begingroup$ Chriserke, you don't want your real-time operation cast your double to int32_t. Scale up (and round) your coefficients (by a factor or 16384) first. In the real-time operation, it should be only int32_t times int32_t into an int32_t result. and that result should be rounded (by adding 8192 or add the bits dropped in the previous sample, if you want this quantization noise-shaped) and shifted right by 14 bits. is your compiler smart enough that the division operation y2 = y2/a10; is actually y2 >>= 14;? $\endgroup$ – robert bristow-johnson May 17 at 19:13
  • $\begingroup$ when i get an hour later today, i will write for you a simple fixed-point Direct Form 1 biquad sample processing loop in C. it only needs about 8 lines of code. $\endgroup$ – robert bristow-johnson May 17 at 19:31
  • $\begingroup$ That would be great. then i will work that into my code and if you would be willing to give it another look over, then i'd be delighted. Kind regards Chris. $\endgroup$ – Chriserke May 17 at 22:23
2
$\begingroup$

A few points to consider

  1. Fixed point DSP is hard. Unless you absolutely cant afford it, try to do whole thing in 32-bit floating point
  2. Assuming you can't you need to manage clipping results and state variable, quantization noise, limit cycles, coefficient quantization, etc.
  3. You do this be carefully choosing representation of in terms number of bits and scaling factor for each individual variable. A good method for this the "Q" notation https://en.wikipedia.org/wiki/Q_(number_format) (which I think RBJ calls)
  4. C is a not a great language for this since its lacking native support for the most basic fixed point multiply function. Most embedded processor do support them through some sort of intrinsic or fake library call. Grab your processor and/or compiler manual and look for things like "fractional multiply", "fractional data" "rounding behavior" "overlfow/underflow" "accumulator" behavior etc.
  5. Understand your requirements: what noise level can you tolerate, how sensitive are too occasional clipping, does the passband actually need to be at 0dB or can you downscale a bit? What's your MIPS budget ? The exact requirements will have significant impact and how you would write the code. Good news is: you really don't need much memory at all.
  6. Your specific example: with 16-bits for the coefficients, you are getting a fair bit of coefficient quantification especially on the zeros. This means that your cutoff frequency and passband gain will be slightly off, which may or may not be a problem
  7. Your code as written isn't really fixed point. You keep the state variables as doubles and so the compiler will simply promote all ints to doubles before the multiply and then truncate the final result down to int_16 when assigning to the output variable.
  8. Get rid of the divides. They don't do anything useful.
  9. "If" statements can be VERY expensive on processors with a deep pipeline. Try to manage clipping through proper scaling of all variables and avoid the "if" statements all together.
  10. First thing I would look for is whether the processor or compiler supports a 32-bit fractional multiply at a responsable cost. This multiplies two 32-bit numbers and keeps 32 MSBs of the 64-bit result. Then I'd keep the coefficients in Q30.1, the state variables in Q28.3 and scale the final result back down to Q15.0.
  11. If available, set rounding mode to "round towards zero", that's slightly more noisy but guarantees to eliminate any limit cycles.

More good reading

https://en.wikipedia.org/wiki/Fixed-point_arithmetic

$\endgroup$
  • $\begingroup$ Hello, This is done on a Xmega chip where 32bits take about 4 times the amount of time a 16 bit takes with no floating point processor. As for passband gain its specified at +1 and -4 throughout the entire passband and -3db at 500hz.As for the states, i should change those to ints as well in the 3rd example code? Lastly a designated dsp is not used but a all in 1 microcontroller with adc and dac. $\endgroup$ – Chriserke May 18 at 15:36
  • $\begingroup$ Does the processor support fractional data types and fractional multiplies? I suggested you keep the states in Q28.3. If you don't know what that is, it's probably best to study fixed point math for a bit before starting to write code. You'll need it. I assume the process has a floating point library. Try writing it all in 32-bit floats and benchmark it. If that fits your MIPS budget, its by far the easiest way to go $\endgroup$ – Hilmar May 18 at 15:43
  • $\begingroup$ The very first code i posted is in floating point which does not reach the speed. Thus going for fixed point in which shifting is faster than multiplying for my microcontroller i believe. Also there are no libraries for math used $\endgroup$ – Chriserke May 18 at 15:52
  • $\begingroup$ I understand what for example Q1.14 is or 28.3 its 1int bit 14 fractions and 28int bits3 fractions i don't know how to correctly implement it however which is why i am struggling $\endgroup$ – Chriserke May 18 at 17:22
  • $\begingroup$ Your first code is in 64-bit double precision. Turn it into 32-bit single precision (use data type "float"), unroll the loops, get rid of the divides and check whether this is fast enough. If not, check whether your processor or compiler support fractional multiplies $\endgroup$ – Hilmar May 18 at 21:31
0
$\begingroup$

If you multiply 2 16-bit numbers (int16_t or S.15), the result requires a 32-bit variable (int32_t or S.31), and the 2 pre-scaling values, if any, add.

$\endgroup$
  • $\begingroup$ Hello, Thanks for the tip. I've added my fixed point attempt. i used 16bits ints and multipled my coefficients by 2^14 because with 2^15 it won't fit in a 16 bits int. Would you mind having a look to see if this is what u meant? King regards Chris $\endgroup$ – Chriserke May 17 at 18:26
  • $\begingroup$ hot, he ain't doing S.31. he's doing S17.14 . also C does not automatically extend word width. we gotta be careful. Chriserke needs to scale up all coefficients by $2^{14}$ and do all of the fixed-point arithmetic in int32_t and scale the result down by $2^{-14}$ by shifting the result to the right by 14 bits. some form of rounding (or better yet, fraction saving or noise shaping) should be done before shifting right by 14 bits. $\endgroup$ – robert bristow-johnson May 17 at 19:17
  • $\begingroup$ @robertbristow-johnson i've sacled up the coefficients in my edit where i included the fixed point attempt. However i don't think i've done the actual arithmic in 32 bit but in 16 bit instead. However im unfamiliar with rounding, fraction saving and noise shaping. $\endgroup$ – Chriserke May 17 at 19:31
  • $\begingroup$ but xA[0] is a double, not an int32_t. what's going on there? $\endgroup$ – robert bristow-johnson May 17 at 19:32
  • $\begingroup$ you know what rounding is. but i will demonstrate a simple fraction saving (which is first-order noise shaping) and with only 16/32 bit precision, you're gonna need some noise shaping if you wanna eliminate some limit cycles and other nasties. $\endgroup$ – robert bristow-johnson May 17 at 19:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.