How to interpret integer coefficients for a filter generated by Matlab

Question

I have always used high-performance 32-bit microcontrollers with hardware FPU support. To calculate the coefficients of digital filters, I used the fdatool package in Matlab. I used single-precision floating point coefficients for digital filters. Next, in accordance with this figure, I implemented digital filters: Moreover, sometimes it was possible to do a little optimization, for example with the LP filter I did this:

void lp_filter (float *x ,float *y,float *aa ,float g ,float vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=vh*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0] + x[1] + x[1] + x[2] - aa[0]*y[1] - aa[1]*y[2]); //b0=1 b1=2 b2=1
}

I knew that the b coefficients would always be (1,2,1), so I simply could not do unnecessary multiplication operations. Perhaps the compiler itself could do this, but I'm not sure about that...

Now I have a low-performance 32-bit microcontroller without an FPU. I wanted to implement digital filters, but with integer mathematics. The fdatool package also allows you to calculate integer filter coefficients, but I don’t understand how to use them... For example, here are the coefficients for a 500Hz 6th order low pass filter:

#define MWSPT_NSEC 7
const int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1 };
const int32_T NUM[MWSPT_NSEC][3] = {
  {
     265839664,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
     209661133,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
     186862344,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
    2147483647,           0,           0 
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1,3,1,3,1 };
const int32_T DEN[MWSPT_NSEC][3] = {
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -2147483648,  1483049895 
  },
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -2024667000,   715827883 
  },
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -1804502424,   404468153 
  },
  {
    2147483647,           0,           0 
  }
};

I choose the format of the filter coefficients as 32 bit signed integer. If I follow the same implementation I used for floating point, the input value may overflow 32 bits... Perhaps I need to use the coefficients as 16 bit signed integers or store the results in 64 bit registers? I tried using the int64_t format, but the filters don't work correctly. I can't understand why... Also, in the integer implementation of digital filters, I don’t understand how I should work with the final result? For example, how can I understand where the comma would be if it were a float number?

---

For digital filters I used single precision floating point coefficients. Here is an example of such coefficients that I obtained in Matlab:

#define MWSPT_NSEC 3
const int NL[MWSPT_NSEC] = { 1,3,1 };
const real32_T NUM[MWSPT_NSEC][3] = {
  {
    0.01185768284,              0,              0 
  },
  {
                1,              2,              1 
  },
  {
                1,              0,              0 
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1 };
const real32_T DEN[MWSPT_NSEC][3] = {
  {
                1,              0,              0 
  },
  {
                1,   -1.669203162,   0.7166338563 
  },
  {
                1,              0,              0 
  }
};

Based on these coefficients, I implemented a digital filter like this:

float sek1_x[3],sek1_y[3],sek1_a[3]={1,   -1.669203162,   0.7166338563},sek1_b[3]={1,              2,              1};
void lp_filter (float *x ,float *y,float *aa ,float *bb ,float g ,float in)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=in*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0]*bb[0] + x[1]*bb[1] + x[2]*bb[2] - aa[1]*y[1] - aa[2]*y[2]);
} 
//Calling the filter function:
lp_filter(sek1_x,sek1_y,sek1_a,sek1_b,0.01185768284,adc_res);

The fdatool package also allows you to calculate integer filter coefficients. Here is an example of such coefficients for exactly the same filter, which I obtained in Matlab:

#define MWSPT_NSEC 3
const int NL[MWSPT_NSEC] = { 1,3,1 };
const int32_T NUM[MWSPT_NSEC][3] = {
  {
      25464180,           0,           0 
  },
  {
    2147483647,  2147483647,  2147483647 
  },
  {
    2147483647,           0,           0 
  }
};
const int DL[MWSPT_NSEC] = { 1,3,1 };
const int32_T DEN[MWSPT_NSEC][3] = {
  {
    2147483647,           0,           0 
  },
  {
    2147483647, -2147483648,  1538959525 
  },
  {
    2147483647,           0,           0 
  }
};

Here is the code for implementing a digital filter with integer coefficients based on my code (as well as the figure I provided above) for implementing floating point filters:

int32_t sek1_x[3], sek1_y[3];
int32_t sek1_a[3] = { 2147483647, -2147483648, 1538959525 };
int32_t sek1_b[3] = { 2147483647, 2147483647, 2147483647 };

void lp_filter(int32_t *x, int32_t *y, int32_t *aa, int32_t *bb, int32_t g, int32_t in) {
    x[2] = x[1];
    x[1] = x[0];
    x[0] = ((int64_t)in * g) >> 31; // Scaling by shifting

    y[2] = y[1];
    y[1] = y[0];
    
    // Perform the filtering with scaling
    y[0] = (((int64_t)x[0] * bb[0] + (int64_t)x[1] * bb[1] + (int64_t)x[2] * bb[2] - (int64_t)aa[1] * y[1] - (int64_t)aa[2] * y[2])*2147483647) >> 31;
}

// Calling the filter function:
lp_filter(sek1_x, sek1_y, sek1_a, sek1_b, 25464180, adc_res);

I'm not sure that this will work, since I don't understand what part of the data is needed for, namely these lines:

const int DL[MWSPT_NSEC] = { 1,3,1 };
const int32_T DEN[MWSPT_NSEC][3] = {
  {
    2147483647,           0,           0 <- this line
  },
  {
    2147483647, -2147483648,  1538959525 
  },
  {
    2147483647,           0,           0 <- this line
  }
};

I also don’t understand how I can properly monitor overflow. I did shifts to try to avoid overflow, but I'm not sure if it's correct... Does anyone have any ideas on how to correctly use integer coefficients generated by matlab to implement digital filters?

---

I wrote an example from the information that I understood, but I made a mistake somewhere, but I don’t understand where... Source code with float filter implementation (6th order low pass filter and 2nd order high pass filter):

float   lp500_sek1_x[3],lp500_sek1_y[3],lp500_sek1_a[2]={-1.833933353,   0.8860449791},
        lp500_sek2_x[3],lp500_sek2_y[3],lp500_sek2_a[2]={-1.669203162,   0.7166338563},
        lp500_sek3_x[3],lp500_sek3_y[3],lp500_sek3_a[2]={-1.586906791,   0.6319990754},
        hp5_sek1_x[3],hp5_sek1_y[3],hp5_sek1_a[2]={-1.998889327,   0.9988899231};

//Calling filters
lp_filter(lp500_sek1_x,lp500_sek1_y,lp500_sek1_a,0.01302789338,((float)res_adc));
lp_filter(lp500_sek2_x,lp500_sek2_y,lp500_sek2_a,0.01185768284,lp500_sek1_y[0]);
lp_filter(lp500_sek3_x,lp500_sek3_y,lp500_sek3_a,0.01127306558,lp500_sek2_y[0]);
hp_filter(hp5_sek1_x,hp5_sek1_y,hp5_sek1_a,0.9994447827,lp500_sek3_y[0]);

//Filter functions
void lp_filter (float *x ,float *y,float *aa ,float g ,float vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=vh*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0] + x[1]+x[1] + x[2] - aa[0]*y[1] - aa[1]*y[2]); //b0=1 b1=2 b2=1
}
void hp_filter (float *x ,float *y,float *aa ,float g ,float vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=vh*g;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=(x[0] - (x[1]+x[1]) + x[2] - aa[0]*y[1] - aa[1]*y[2]); //b0=1 b1=-2 b2=1
}

Source code with an integer implementation of the same filters (6th order low-pass filter and 2nd order high-pass filter):

int32_t ilp500_sek1_a[2]={-1969170942, 951383551},
        ilp500_sek2_a[2]={-1792293247, 769479743},
        ilp500_sek3_a[2]={-1703928191, 678603839},
        ihp5_sek1_a[2]={-2146291070,  1072549887},
        ilp6_sek1_a[2]={-2145098367, 1071359231};

int64_t ilp500_sek1_x[3],ilp500_sek1_y[3],
        ilp500_sek2_x[3],ilp500_sek2_y[3],
        ilp500_sek3_x[3],ilp500_sek3_y[3],
        ihp5_sek1_x[3],ihp5_sek1_y[3],
        ilp6_sek1_x[3],ilp6_sek1_y[3];

//Calling filters
lp_filter_int(ilp500_sek1_x,ilp500_sek1_y,ilp500_sek1_a,13988593,((int64_t)res_adc));
lp_filter_int(ilp500_sek2_x,ilp500_sek2_y,ilp500_sek2_a,12732089,ilp500_sek1_y[0]);
lp_filter_int(ilp500_sek3_x,ilp500_sek3_y,ilp500_sek3_a,12104361,ilp500_sek2_y[0]);
hp_filter_int(ihp5_sek1_x,ihp5_sek1_y,ihp5_sek1_a,1073145662,ilp500_sek3_y[0]);

//Filter functions
void lp_filter_int (int64_t *x ,int64_t *y,int32_t *aa ,int32_t g ,int64_t vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=(vh*(int64_t)g)>>30;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=((x[0]*(int64_t)0x3FFFFFFF) + (x[1]*(int64_t)0x7FFFFFFF) + (x[2]*(int64_t)0x3FFFFFFF) - ((int64_t)aa[0]*y[1]) - ((int64_t)aa[1]*y[2]))>>30; //b0=1 b1=2 b2=1
}
void hp_filter_int (int64_t *x ,int64_t *y,int32_t *aa ,int32_t g ,int64_t vh)
{
    x[2]=x[1];
    x[1]=x[0];
    x[0]=(vh*(int64_t)g)>>30;
    y[2]=y[1];
    y[1]=y[0];
    y[0]=((x[0]*(int64_t)0x3FFFFFFF) - (x[1]*(int64_t)0x7FFFFFFF) + (x[2]*(int64_t)0x3FFFFFFF) - ((int64_t)aa[0]*y[1]) - ((int64_t)aa[1]*y[2]))>>30; //b0=1 b1=-2 b2=1
}

After the integer low-pass filter, I see a similar signal as after the float low-pass filter, but after the high-pass filter everything breaks down... I can’t understand what I’m doing wrong?

Answer 1

There is a lotta detail in the question.

I have done a hella lotta fixed-point DSP in my day. Both with a fixed-point DSP (56K) and with an integer CPU such as the 68K but also in generic C code. I posted an answer about specific C code that does a biquad in fixed point and uses noise-shaping at the sole quantization point. You should look at it.

Here is the main principle. You saw a comment refer to the "Q" notation. E.g. Q2.14 means a 16-bit fixed-point number with 2 bits to the left of the binary point and 14 to the right. What this really means is that the 16-bit integer has an integer value (this is the value that the C program sees) that is scaled up by a factor of $2^{14}$. This doesn't really matter for the signal; whatever scaling of the signal in the input is the same scaling for the output. But it does matter for the coefficients. Every coefficient is represented as an integer that is $2^{14}$ times larger than the actual numerical value of the coefficient.

So, when you integer multiply that coefficient times another value (like a signal or a state), the result will be $2^{14}$ times bigger than it should be. So you'll have to shift it 14 bits to the right to scale it with $2^{-14}$ and undo the scaling of the coefficient necessary to make it an integer. This shifting right causes 14 bits to fall offa the edge, which is quantization error. Noise shaping can help minimize the consequences of that quantization.

With 32 bits, the quantization error will be much smaller in effect.

So you have to deal with the range of your coefficients, so you know how much to scale them. For a biquad, we know the denominator coefficients will have absolute value of less than 2. So we can represent those coefficients well as Q2.30 fixed-point values. So using int32_t as your coefficients and int64_t as your accumulator, then after multiplying, you will have to shift the 64-bit value 30 bits to the right.

Besides noise-shaping (in this specific case, "fraction saving" is the noise shaping), there are some other useful tricks in fixed-point DSP. Another trick is to use the Direct Form 1 and do not use Direct Form II. The reasons why are two-fold: First the DF1 has only one point of quantization. And second, the DF1 puts the effect of the zeros before the poles, so in a high-Q filter, there is attenuation (from a zero) that will occur before there is gain (from the pole). This means that unnecessary internal overflow is avoided. And keeping a 64-bit accumulator and adding up all 5 values before quantization allows you to quantize after the addition.