Floating point to fixed point in C

Question

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);

Answer 1

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

Answer 2

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.